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An overview of some analytical approaches to the computation of the struc- 
tural and thermodynamic properties of single component and multicompo- 
nent hard-sphere fluids is provided. For the structural properties, they yield 
a thermodynamically consistent formulation, thus improving and extending 
the known analytical results of the Percus-Yevick theory. Approximate ex- 
pressions for the contact values of the radial distribution functions and the 
corresponding analytical equations of state are also discussed. Extensions of 
this methodology to related systems, such as sticky hard spheres and square- 
well fluids, as well as its use in connection with the perturbation theory of 
fluids are briefly addressed. 

1 Introduction 

In the statistical thermodynamic approach to the theory of simple liquids, 
there is a close connection between the thermodynamic and structural prop- 
erties [1-4]. These properties depend on the intermolecular potential of the 
system, which is generally assumed to be well represented by pair interactions. 
The simplest model pair potential is that of a hard-core fluid (rods, disks, 
spheres, hyperspheres) in which attractive forces are completely neglected. In 
fact, it is a model that has been most studied and has rendered some analytical 
results, although up to this day no general (exact) explicit expression for the 
equation of state is available, except for the one-dimensional case. Something 
similar applies to the structural properties. An interesting feature concerning 
the thermodynamic properties is that in hard-core systems the equation of 
state depends only on the contact values of the radial distribution functions. 
In the absence of a completely analytical approach, the most popular methods 
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to deal with both kinds of properties of these systems are integral equation 
theories and computer simulations. 

It is well known that in real gases and liquids at high temperatures the 
state and thermodynamic properties are determined almost entirely by the 
repulsive forces among molecules. At lower temperatures, attractive forces 
become significant, but even in this case they affect very little the configu- 
ration of the system at moderate and high densities. These facts are taken 
into account in the application of the perturbation theory of fluids, where 
hard-core fluids are used as the reference systems in the computation of the 
thermodynamic and structural properties of real fluids. However, successful 
results using perturbation theory arc rather limited due to the fact that, as 
mentioned above, there are in general no exact (analytical) expressions for 
the thermodynamic and structural properties of the reference systems which 
are in principle required in the calculations. On the other hand, in the realm 
of soft condensed matter the use of the hard-sphere model in connection, for 
instance, with sterically stabilized colloidal systems is quite common. This is 
due to the fact that nowadays it is possible to prepare (almost) monodispersc 
spherical colloidal particles with short-ranged harshly repulsive interparticle 
forces that may be well described theoretically with the hard-sphere potential. 

This chapter presents an overview of the efforts wc have made over the 
last few years to compute the thermodynamic and structural properties of 
hard-core systems using relatively simple (approximate) analytical methods. 
It is structured as follows. In Section 2 we describe our proposals to derive 
the contact values of the radial distribution functions of a multicomponent 
mixture (with an arbitrary size distribution, either discrete or continuous) of 
d-dimensional hard spheres from the use of some consistency conditions and 
the knowledge of the contact value of the radial distribution function of the 
corresponding single component system. In turn, these contact values lead to 
equations of state both for additive and non-additive hard spheres. Some con- 
sequences of such equations of state, in particular the demixing transition, are 
briefly analyzed. This is followed in Section 3 by the description of the Ratio- 
nal Function Approximation method to obtain analytical expressions for the 
structural quantities of three-dimensional single component and multicompo- 
nent fluids. The only required inputs in this approach are the contact values of 
the radial distribution functions and so the connection with the work of the 
previous section follows naturally. Structural properties of related systems, 
like sticky hard spheres or square-well fluids, that may also be tackled with 
the same philosophy are also discussed in Section 4. Section 5 provides an 
account of the reformulation of the perturbation theory of liquids using the 
results of the Rational Function Approximation method for a single compo- 
nent hard-sphere fluid and its illustration in the case of the Lennard-Joncs 
fluid. In the final section, we provide some perspectives of the achievements 
obtained so far and of the challenges that remain ahead. 
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2 Contact Values and Equations of State for Mixtures 

As stated in the Introduction, a nice feature of hard-core fluids is that the 
expressions of all their thermodynamic properties in terms of the radial distri- 
bution functions (RDF) arc particularly simple. In fact, for these systems the 
internal energy reduces to that of the ideal gas and in the pressure equation 
it is only the contact values rather than the full RDF which appear explicitly. 
In this section wc present our approach to the derivation of the contact values 
of hard-core fluid mixtures in d dimensions. 

2.1 Additive Systems in d Dimensions 

If Gij denotes the distance of separation at contact between the centers of 
two interacting fluid particles, one of species i and the other of species j, 
the mixture is said to be additive if cTy is just the arithmetic mean of the 
hard-core diameters of each species. Otherwise, the system is non-additive. 
We deal in this subsection and in Subsection 2.2 with additive systems, while 
non-additive hard-core mixtures will be treated in Subsection 2.3. 

Definitions 

Let us consider an additive mixture of hard spheres (HS) in d dimensions with 
an arbitrary number N of components. In fact, our discussion will remain valid 
for N ^ oo, i.e., for polydisperse mixtures with a continuous distribution of 
sizes. 

The additive hard core of the interaction between a sphere of species i and 

a sphere of species j is (T,y = ■^{(Ji + (Tj), where the diameter of a sphere of 
species i is = . Let the number density of the mixture be p and the mole 
fraction of species the Xi= Pi/ p, where pi is the number density of species i. 
From these quantities one can define the packing fraction 77 = VdpMd, where 
Vd = (77/4) ''/^/r'(l -t- d/2) is the volume of a d-dimensional sphere of unit 
diameter and 

N 

M„^(a")=^x,ar (1) 

i=l 

denotes the nth moment of the diameter distribution. 

In a HS mixture, the knowledge of the contact values gij{(Jij) of the RDF 
gij{r), where r is the distance, is important for a number of reasons. For 
example, the availability of gij{aij) is sufficient to get the equation of state 
(EOS) of the mixture via the virial expression 



2d-i N 



(2) 
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where Z = p/ pksT is the compressibiUty factor of the mixture, p being the 
pressure, fc^ the Boltzmann constant, and T the absohitc temperature. 

The exact form of gij{crij) as functions of the packing fraction 77, the set 
of diameters {ufc}, and the set of mole fractions {xk} is only known in the 
one-dimensional case, where one simply has [5] 

9^J{<^^j) ^ (C^^l)- (3) 

Consequently, for d > 2 one has to resort to approximate theories or empirical 
expressions. For hard-disk mixtures, an accurate expression is provided by 
Jenkins and Mancini's (JM) approximation [6,7], 

% K) = 13^ + 1^ (13^^^^' (^=2)- (4) 
The associated compressibility factor is 

In the case of three-dimensional systems, some important analytical expres- 
sions for the contact values and the corresponding compressibility factor also 
exist. For instance, the expressions which follow from the solution of the 
Percus-Yevick (PY) equation of additive HS mixtures by Lebowitz [8] are 



PY/ ^ 1,3?? <Jt(yjM2 



^ , , 1 M1M2 37? Ml 37?2 , 



Also analytical are the results obtained from the Scaled Particle Theory (SPT) 
[9-12], 



„SPT/_ ^ _ ,3 7? aiajM2 , 3 /^ ajajMa ^ 

9ii i'^^j) - 1 _ ^ + 2 (1 - 7?)2 a,,Ms + 4 (1 - r?)3 1, a,,Ms j ' " 

(8) 

^ / X 1 M1M2 377 M| 37?2 , , 



1 M1M2 37/ M| 37?2 
1^ ^ M3 (1 - 77)2 ^ M| (1 - 77)3 ' 

Neither the PY nor the SPT lead to particularly accurate values and so 
Boublfk [13] and, independently, Grundke and Henderson [14] and Lee and 
Levesque [15] proposed an interpolation between the PY and the SPT contact 
values, that we will refer to as the BGHLL values: 

„BGHLL/, ^_ 1 ,3 V ai(JjM2 rf f aiajM2 Y / , _ 

^""''^ - 1 - 7? + 2 (1 - 77)2 a.jMs +2 (1 - 77)3 1, a,,M, j ' " "^^ 

(10) 
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This leads through Eq. (2) to the widely used and rather accurate Boublik- 
Mansoori-Carnahan-Starhng-Leland (BMCSL) EOS [13,16] for HS mixtures: 

1 , M1M2 3v , M| r?^(3-r?) 
^BMCSL(r?) = + ^^_^^^ +—^^—^, id = 3). (11) 

Refinements of the BGHLL values have been subsequently introduced, among 
others, by Henderson et al. [17-22], Matyushov and Ladanyi [23], and Barrio 
and Solana [24] to eliminate some drawbacks of the BMCSL EOS in the so- 
called colloidal limit of binary HS mixtures. On a different path, but also 
having to do with the colloidal limit, Viduna and Smith [25] have proposed 
a method to obtain contact values of the RDF of HS mixtures from a given 
EOS. However, none of these proposals may be easily generalized so as to be 
valid for any dimensionality and any number of components. Therefore, if one 
wants to have a more general framework able to deal with arbitrary d and N 
an alternative strategy is called for. 



Universality Ansatz 

In order to follow our alternative strategy, it is useful to make use of exact 
limit results that can help one in the construction of approximate expressions 
for gij{<7ij). Let us consider first the limit in which one of the species, say «, 
is made of point particles, i.e., ctj 0. In that case, gii{(Ji) takes the ideal 
gas value, except that one has to take into account that the available volume 
fraction is 1 — 77. Thus, 

\\Tagii{ai) = - -. (12) 

An even simpler situation occurs when all the species have the same size, 
{<^k} o", so that the system becomes equivalent to a single component 
system. Therefore, 

lim gij{(Jij) = gs, (13) 

where gs is the contact value of the RDF of the single component fluid at 
the same packing fraction r] as that of the mixture. Table 1 lists some of 
the most widely used proposals for the contact value gs and the associated 
compressibility factor 

Zs = l + 2'^-^r,gs (14) 

in the case of the single component HS fluid. 

Equations (12) and (13) represent the simplest and most basic conditions 
that gij (aij ) must satisfy. There is a number of other less trivial consistency 
conditions [11, 17, 19, 20, 23, 24, 32 34], some of which will be used later on. 

In order to proceed, in line with a property shared by earlier proposals 
[see, in particular, Eqs. (4), (6), (8), and (10)], we assume that, at a given 
packing fraction r], the dependence of gij{<Tij) on the parameters {cfc} and 
{xfc} takes place only through the scaled quantity 
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Table 1. Some expressions of and Zs for the single component HS fluid. In the 

SHY proposal, r/cp ~ (\/3/6)7r is the crystaUino close-packing fraction for hard disks. 
In the LM proposal, 63 and 64 are the (reduced) third and fourth virial coefficients, 
C(??) = 1.2973(59) -0.062(13)r;/j7cp for d = 4, and <(??) = 1.074(16) + 0.163(45)??/77cp 
for d = 5, where the values of the close-packing fractions are rjcp = ~ 0.617 

and 77cp = 7r^-\/2/30 ~ 0.465 for c? = 4 and d = 5, respectively. 







Zb 


Label Ref. 




1 — (T]/ ID 
(1 - r;)2 


i + ?7 

(1-77)2 


H 


[26] 




1 - r?(277cp - 1)/27j!p 


1 


SHY 


[27] 


1 - 


-2r7 + r?2(,,ep-l)/2r?|p 


l- 277 + 772(77cp-l)/277|p 




27(1 -,7)4 


" 26(1-77)4 


L 


[28] 




1 + 77/2 

(1 - r;)2 


1 + 277 + 3772 

(1 - 77)2 


PY 


[29] 




1-77/2 + 7774 

(1-77)3 


1 + 77 + 77^ 

(1 - 77)3 


SPT 


[9] 




1-77/2 

(1 - r;)3 


1 + 77 + 77^-77^ 


CS 


[30] 


1 


f [2'-%, - av)b4/b-s]ri 


l+2^'-S,flLM 


LM 


[31] 


i - a'/) 


(''4/6:i)'/-[C('/)-l]2-^ ''/-4/r 



CTij Aid 

More specifically, we assume 

9ij{(rij) =Qiv,Zij), (16) 

where the function ^(77, z) is universal in the sense that it is a common func- 
tion for all the pairs {i, j) , regardless of the composition and number of compo- 
nents of the mixture. Of course, the function Q{r], z) is in principle different for 
each dimensionality d. To clarify the implications of this universality ansatz, 
let us imagine two mixtures Ad and M' having the same packing fraction 
77 but strongly differing in the set of mole fractions, the sizes of the parti- 
cles, and even the number of components. Suppose now that there exists a 
pair in mixture M. and another pair in mixture A4' such that 

Zij = Zi'j' . Then, according to Eq. (16), the contact value of the RDF for the 
pair in mixture M. is the same as that for the pair in mixture 

A4' , i.e., gij{(7ij) = (ji'j'{<yi'j')- In order to ascribe a physical meaning to the 
parameter Zij, note that the ratio M^-i/M^ can be understood as a "typ- 
ical" inverse diameter (or curvature) of the particles of the mixture. Thus, 
z~^ = + cr~^)/{Mci^i/Mfi) represents the arithmetic mean curvature, 

in units of M^-i/Md, of a particle of species i and a particle of species j. 
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Once the ansatz (16) is adopted, one may use the limits in (12) and (13) 
to get Q{rj,z) at z = and z = 1, respectively. Since Zn ^ in the limit 
fTj 0, insertion of Eq. (12) into (16) yields 

g{^^Q) = J—^g^{r,). (17) 

1 - ri 

Next, if all the diameters are equal, Zij 1, so that Eq. (13) implies that 

Giv, 1) = 5s. (18) 



Linear Approximation 

As the simplest approximation [35], one may assume a linear dependence of 
g on z that satisfies the basic requirements (17) and (18), namely 

Inserting this into Eq. (16), one has 

el/ \ ^ , f ^ \ Md-l (Tiaj 

4K-) = T3^+(5-T3^j^^- (20) 

Here, the label "el" is meant to indicate that (i) the contact values used are 
an extension of the single component contact value Qs and that (ii) Q{r},z) is 
a linear polynomial in z. This notation will become handy below. Although 
the proposal (20) is rather crude and does not produce especially accurate 
results for gij{aij) when d > 3, it nevertheless leads to an EOS that exhibits 
an excellent agreement with simulations in 2, 3, 4, and 5 dimensions, provided 
that an accurate (?s is used as input [35-39]. This EOS may be written as 

Zei(77) = 1 + -^2^-^00 - Q,) + [Z,{r]) - 1] Ou (21) 
l-rj 

where the coefficients i7m depend only on the composition of the mixture and 
are defined by 

= 2-^'^-)^ E ('-7)m„+„M,_„. (22) 

d n=0 ^ ' 

In particular, for d = 2 and d = 3, 



1 Ml 



1-TJ M2 



1 



{d = 2), (23) 
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Zeliv) 



- V 
3?7 



1 



M1M2 
2M3 

1 - 



Ml 



M1M3 



1-7/ 

(d = 3). 



1 + 



M1M3 



(24) 



As an extra asset, from Eq. (21) one may write the virial coefficients of 
the mixture B„, defined by 



(25) 



in terms of the (reduced) virial coefficients of the single component fluid &„ 
defined by 



The result is 



(26) 



(27) 



In the case of binary mixtures, these coefficients are in very good agreement 
with the available exact and simulation results [35,37], except when the mix- 
ture involves components of very disparate sizes, especially for high dimen- 
sionalities. One may perform a slight modification such that this deficiency is 
avoided and thus get a modified EOS [37, 40] . For d = 2 and d = 3 it reads 



1 

+X2 



1 - 



m 



1 - r?i 



1 - »?2 

- 



Zsiv) 



(72 — CTi 



<72 



d-1 



<7l — 02 



(d = 2,3). 



(28) 



where r\i = vapicrf is the partial volume packing fraction due to species i. In 
contrast to most of the approaches (PY, SPT, BMCSL, el, . . . ), the proposal 

(28) expresses Z{ri) in terms not only of Zs{ri) but also involves Z^ ^3:^^ 

and Zs ^Y^^^. Equation (28) should in principle be useful in particular for 

binary mixtures involving components of very disparate sizes. However, it is 
slightly less accurate than the one given in Eq. (21) for ordinary mixtures [37]. 



Quadratic Approximation 

In order to improve the proposal contained in Eq. (20), in addition to the 

consistency requirements (12) and (13), one may consider the condition stem- 
ming from a binary mixture in which one of the species (say i = 1) is much 
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larger than the other one (i.e., (J\/(J2 oo), but occupies a negligible volume 
(i.e., xi{(Ti/a2y 0). In that case, a sphere of species 1 is felt as a wall by 
particles of species 2, so that [17,20,41] 



lim [9i2{(Ti2) - 2"^ ^77522(0-2)] = 1. 



(29) 



Hence, in the limit considered in Eq. (29), we have Z22 ^ 1, -212 ^ 2. Conse- 
quently, under the universality ansatz (16), one may rewrite Eq. (29) as 



g(77,2) = l + 2'^-Se(?7,l) 



(30) 



Thus, Eqs. (17), (18), and (30) provide complete information on the function 
G at z = 0, z = 1, and z = 2, respectively, in terms of the contact value Qs of 
the single component RDF. 

The simplest functional form of G that complies with the above consistency 
conditions is a quadratic function of z [42]: 



Giv, z) = Goiv) + Giiv)z + g2{v)z^ 

where the coefficients Giirj) and G2{v) ^re exphcitly given by 

2 -,7/2 



V 



l-V/2 



1 



V 



(1 



-,d-2 



V)9s- 



(31) 

(32) 

(33) 



Therefore, the explicit expression for the contact values is 



1 



1-7J 
'1 - 



+ 



(2 
77/2 



id-2 



V)9s 



2-77/2' 



Md (Tij 



1-77 



- (1 - 2''-^n)9, 



Md 



(34) 



Following the same criterion as the one used in connection with Eq. (20), the 

label "e2" is meant to indicate that (i) the resulting contact values represent 
an extension of the single component contact value Qs and that (ii) G{ri^z) 
is a quadratic polynomial in z. Of course, the quadratic form (31) is not the 
only choice compatible with conditions (17), (18), and (30). For instance, a 
rational function was also considered in Ref. [42]. However, although it is 
rather accurate, it does not lead to a closed form for the EOS. In contrast, 
when Eq. (34) is inserted into Eq. (2), one gets a closed expression for the 
compressibility factor in terms of the packing fraction 77 and the first few 
moments M„, n < d. The result is 



Ze2{v) = 1 + 2 



d-2 



1-77 

+ [ZM - 1] [2f2i 



[2(l2o - 2121 + i^2 



no 



-^d-2 



{^2 



[Qi - U2H 
-^1)77] , 



(35) 
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where the quantities O^n are defined in Eq. (22). Quite interestingly, in the 
two-dimensional case Eq. (35) reduces to Eq. (23), i.e.. 



Z^iirj) = Z^2iv), {d=2). 



(36) 



This illustrates the fact that two different proposals for the contact values 
gij{<Tij) can yield the same EOS when inserted into Eq. (2). On the other 
hand, for three-dimensional mixtures Eq. (35) becomes 



+ 



Ml Ma 



1 



r?- 



1-V ' M3 
which differs from Eq. (24). In fact, 

M1M2 



M1M3 



Zsiv) 



Zeliv) - Ze2{v) 



2M3 



1 



Ml Ms 



1-V 



(1 - 2jj)Z,{r,) 



{d = 3), 
(37) 

{d = 3). 
(38) 



Specific Examples 

In this subsection, rather than carrying out an exhaustive comparison with 
the wealth of results available in the literature, we will consider only a few 
representative examples. In particular, for d = 3, we will restrict ourselves 
to a comparison with classical proposals (say BGHLL, PY, and SPT for the 
contact values). The comparison with more recent ones may be found in Refs. 
[35,42,43]. 

Thus far the development has been rather general since remains free in 
Eqs. (20) and (34). In order to get specific results, it is necessary to fix [cf. 
Table 1]. In the one-dimensional case, one has gs = 1/(1 — 7?) and so one gets 
the exact result (3) after substitution into Eq. (20). Similarly Eqs. (32) and 
(33) lead to Qi = Q2 = ^ so we recover again the exact result. 

If in the two-dimensional case we take Henderson's value [26] gs = gf, 
then the linear approximation (20) reduces to the JM approximation, Eq. (4) . 
This equivalence can be symbolically represented as g°j^^ = gfj^, where the 
label "eHl" refers to the extension of Henderson's single component value in 
the linear approximation. While 5™ is very accurate, even better results are 
provided by the quadratic form (34), especially if Luding's value [28] gs = gl' 
is used [44]. 

In the three-dimensional case, Eq. (20) is of the form of the solution of the 
PY equation [8]. In fact, insertion of gs = gf^ leads to Eq. (6), i.e., 5|_f^^ = 
gfj^- Similarly, if the SPT expression [9] gs = g^"^ is used for the single 
component contact value in the quadratic approximation (34), we reobtain 
the SPT expression for the mixture, Eq. (8). In other words, gf^^^^ = g^f^ ■ 
On the other hand, if the much more accurate CS [30] expression gs = g^^ is 
used as input, we arrive at the following expression: 
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9ir 



1-7/2 (l-r?)2 



7?/3) aiajM2 _^ rf 



(l-r//2) fa^ajA'h 



(1 - 77)3 



(TijMs 



(d = 3), 



which is different from the BGHLL one, Eq. 



(39) 

(10), improves the latter for 



Zij > 1, and leads to similar results for Zij < 1, as comparison with computer 
simulations shows [42]. The four approximations (6), (8), (10), and (39) are 
consistent with conditions (12) and (13), but only the SPT and eCS2 are also 
consistent with condition (29). It should also be noted that if one considers 
a binary mixture in the infinite solute dilution limit, namely Xi 0, so that 
Zi2 ^ 2/(1 + c72/(Ji), Eq. (39) yields the same result for (712(0-12) as the one 
proposed by Matyushov and Ladanyi [23] for this quantity on the basis of 
exact geometrical relations. However, the extension that the same authors 
propose when there is a non- vanishing solute concentration, i.e., for xi ^ 0, 
is different from Eq. (39). 

Equation (34) can also be used in the case of hyperspheres {d > 4) [42]. In 
particular, a very good agreement with available computer simulations [38] is 
obtained for d = 4 and d = 5 by using Luban and Michels [31] value Qs = gl'^- 



N 
IN 



0.08 
0.06 
0.04 
0.02 
0.00 
-0.02 



' 1 ' 

PY 


1 ' 1 . ' 1 ' 1 
/ 

/' ; 


. SPT 


' 


eCS1 


/ d 


eCS2 


/ 




/' c> /• 

. . - ■ ' 


— ^^^■G=-^- 




1 


1 ^ 



0.0 0.1 0.2 0.3 0.4 0.5 



Fig. 1. Deviation of the compressibility factor from the BMCSL value, as a function 
of the packing fraction 77 for an equimolar three-dimensional binary mixture with 
o'2/o'i = 0.6. The open (Ref. [18]) and closed (Rcf [45]) circles are simulation data. 

The lines are the PY EOS ( ), the SPT EOS ( ), the eCSl EOS (•••), and 

the eCS2 EOS ( ). 
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Fig. 2. Compressibility factor for three equimolar mixtures in 4D and 5D systems. 
Lines are the eLMl predictions, while symbols are simulation data [38]. 



Now we turn to the compressibility factors (21) and (35), which are ob- 
tained from the contact values (20) and (34), respectively. Since they depend 
on the details of the composition through the d first moments, they are mean- 
ingful even for continuous poly disperse mixtures. 

As said above, in the two-dimensional case both Eqs. (21) and (35) reduce 
to Eq. (23), which yield very accurate results when a good Zs is used as 
input [39,42,44]. For three-dimensional mixtures, insertion of Zs = Z^^ in 
Eqs. (24) and (37) yields 

^eCSl(r?)=^BMCSL(j?) + ^^-f^P^^(MiM3-M2), (d = 3), (40) 
^eCS2(»y) = ^BMCSL(7?)-^^-3^P^^(MiM3-M|), (d = 3), (41) 

where Zbmcsl(''/) is given by Eq. (11). Note that Zccsi(»7) > ^bmcsl(?7) > 
■^eCS2(??)- Since simulation data indicate that the BMCSL EOS tends to un- 
derestimate the compressibility factor, it turns out that, as illustrated in Fig. 
1 for an equimolar binary mixture with a^j cr\ — 0.6, the performance of ^ccsi 
is, paradoxically, better than that of .Z'ecs2 [42] , despite the fact that the un- 
derlying linear approximation for the contact values is much less accurate than 
the quadratic approximation. This shows that a rather crude approximation 
such as Eq. (20) may lead to an extremely good EOS [35,37-39], which, as 
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clearly seen in Fig. 1, represents a substantial improvement over the classical 
proposals. Interestingly, the EOS corresponding to ^ocsi has recently been 
independently derived as the second order approximation of the Fundamental 
Measure Theory for the HS fluid by Hansen-Goos and Roth [46]. 

In the case of = 4 and d = 5, use of Zsirj) = ^s^'^(^) in Eq. (21) 
produces a simple extended EOS of a mixture of hard additive hyperspheres 
in these dimensionalities. The accuracy of these two EOS for hard hypersphere 
mixtures in the fluid region has been confirmed by simulation data [38] for a 
wide range of compositions and size ratios. In Fig. 2, this accuracy is explicitly 
exhibited in the case of three equimolar mixtures, two in 4D and one in 5D. 

2.2 A More Consistent Approximation for Three-Dimensional 
Additive Mixtures 

Up to this point, we have considered an arbitrary dimensionality d and have 
constructed, under the universality assumption (16), the acurate quadratic 
approximation (34), which fulfills the consistency conditions (12), (13), and 
(29). However, there exist extra consistency conditions that are not necessarily 
satisfied by (34). In particular, when the mixture is in contact with a hard wall, 
the state of equilibrium imposes that the pressure evaluated near the wall by 
considering the impacts with the wall must be the same as the pressure in the 
bulk evaluated from the particle-particle collisions. This consistency condition 
is especially important if one is interested in deriving accurate expressions for 
the contact values of the particle- wall correlation functions. 

Since a hard wall can be seen as a sphere of infinite diameter, the contact 
value g^,j of the correlation function of a sphere of diameter aj with the wall 
can be obtained from gij{aij) as 

gu,j = lim gij{aij). (42) 

17^ — KX> 

XiCrf->0 

Note that g^j provides the ratio between the density of particles of species 

j adjacent to the wall and the density of those particles far away from the 
wall. The sum rule connecting the pressure of the fluid and the above contact 
values is [47] 

N 

Zw{il) = ^Xjgwj, (43) 

where the subscript been used to emphasize that Eq. (43) repre- 

sents a route alternative to the virial one, Eq. (2), to get the EOS of the HS 
mixture. The condition Z = Zw is equivalent to (29) in the special case where 
one has a single fluid in the presence of the wall. However, in the general case 
of a mixture plus a wall, the condition Z = Z^ is stronger than Eq. (29). 
In the two-dimensional case, it turns out that the quadratic approximation 
(34) already satisfies the requirement Z = Z^ , regardless of the density and 
composition of the mixture [44]. However, this is not the case for d > 3. 
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Our problem now consists of computing Qij {aij ) and the associated g^jj for 
the HS mixture in the presence of a hard wall, so that the condition Z = 
is satisfied for an arbitrary mixture [43]. Due to the mathematical complexity 
of the problem, here we will restrict ourselves to three-dimensional systems 
(d — 3). Similarly to what we did in the preceding subsection, we consider 
a class of approximations of the universal type (16), so that conditions (12) 
and (13) lead again to Eqs. (17) and (18), respectively. Notice that Eq. (16) 
implies in particular that 

M2 

9wj = G{v,Zyjj), Zyjj = (44) 

Assuming that z = is a regular point and taking into account condition 
(17), t/(?7, z) can be expanded in a power series in z: 

00 

71=1 

After simple algebra, using the ansatz (16) and Eq. (45) in Eqs. (2) (with 
d= 3) and (43) one gets 

Ml Mo -A Mo" , 

^ = ^° + ^"^^^^ + E E ^^^.^r^"4-". (46) 

n=l 3 iij=l 



Notice that if the series (45) is truncated after a given order n>3,Zyj is given 
by the first n moments of the size distribution only. On the other hand, Z still 
involves an infinite number of moments if the truncation is made after n > 4 
due to the presence of terms like J2i,j ^i^j'^i'^i^ l^iji T^i j ^^i^'j^^^j^l^'ij- ■ ■ ■ ■ 
Therefore, if we want the consistency condition Z = Z^to be satisfied for any 
discrete or continuous polydispcrse mixture, either the whole infinite series 
(45) needs to be considered or it must be truncated after n — 3. The latter is 
of course the simplest possibility and thus we make the approximation 

Qin, z) = Goiv) + Qi{ri)z + g2{v)z^ + Qz{'n)z\ (48) 

As a consequence, Z and Z^ depend functionally on the size distribution of 
the mixture only through the first three moments (which is in the spirit of 
Rosenfeld's Fundamental Measure Theory [48]). 

Using the approximation (48) in Eqs. (46) and (47) we are led to 



Z = gQ + r] 



M1M2 , , Mn , 

(3^0 + 2^1) + 2^ (Si + 202 + 2^3) 



(49) 
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Z^ = g, + 2^^ei + + 2^3) . (50) 

Thus far, the dependence of both Z and Z^, on the moments Mi, M2, and M3 
is explicit and we only lack the packing- fraction dependence of Qi, G2, and 
Q3. From Eqs. (49) and (50) it follows that the difference between Z and Zyj 
is given by 



M1M2 , , . , M? 



= ^nrZ^ [3^^o - 2(1 - V)gi] + 2'i^ [vGi - 2(1 - v)G2 - 2(2 - 7^)03] . 



Therefore, Z = Zyj for any dispersity provided that 



(51) 



S.W = ^, (52) 

where use has been made of the definition of Qq, Eq. (17). To close the problem, 
we use the equal size limit given in Eq. (18), which yields G0+G1+G2+G3 = ^s- 
After a little algebra we are led to 

g,{rj) = {2-v)9s-'^^^^, (54) 

(1 - V) 

g3iv)^{l-ri){9f^~9s). (55) 

This completes the derivation of our improved approximation, which we will 
call "e3" , following the same criterion as the one used to call "el" and "e2" to 
the approximations (20) and (34), respectively. In Eq. (55), gf^""" is the SPT 
contact value for a single fluid, whose expression appears in Table 1. From 
Eq. (55) it is obvious that the choice Qs = gf^'^ makes our eS approximation 
to become the e2 approximation, both reducing to the SPT for mixtures, Eq. 
(8). This means that the SPT is fully internally consistent with the require- 
ment Z = Z^, although it has the shortcoming of not being too accurate in 
the single component case. The c3 proposal, on the other hand, satisfies the 
condition Z = Zyj and has the flexibility of accommodating any desired gg- 

For the sake of concreteness, let us write explicitly the contact values in 
the e3 aproximation: 



2 



^-V 2(1-77)^^3 a., 



2 + r?V4 
(2 - l)9s - 72- 
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e3 1 377 M2 ^ 



1 - ?? (1 - ??) ^3 



(2 - r?)5s - 



2 + 7?V4 



M2 
M3' 



+8(1 - n) {g.. 



SPT 



M2 



(57) 



With the above results the compressibihty factor may be finally written in 
terms of Zg as 



ZM = 



(i-r?) 



M1M2 



Ml 



37? 



1 



(58) 



A few comments are in order at this stage. First, from Eq. (49) we can 
observe that, for the class of approximations (48), the compressibility factor 
Z does not depend on the individual values of the coefBcients Q2 and ^3, 
but only on their sum. As a consequence, two different approximations of 
the form (48) sharing the same density dependence of Qi and Q2 + also 
share the same virial EOS. For instance, if one makes the choice = fff^i 
then Zf,pYz = ZpYi even though fjlj^^icrij) 7^ sfj^i^ij)- Furthermore, if one 
makes the more accurate choice gs — g^^, then Zecsi = .^bmcsl, but again 



The eCS3 contact values are 



+ 



37] 



+ 



V 2{l-r]) Ms (7ij 

,2 __X3 



M2 OiGj ^ rf{l 



4(l-7?)3 VMs Gij 



M2 (JjCTj 

4(l-r?)2 VMs an 



(59) 



,CS3 1 , 377 M2 , r/^(l + 7?) (M2 



t 

(1-77)2 VM3 



r^-. 1 . (60) 



In Figs. 3 and 4 we display the performance of the contact values as given 
by Eqs. (59) and (60), respectively, by comparison with results of computer 
simulations for both discrete and poly disperse mixtures. In both figures we 
have also included the results that follow from the classical proposals as well 
as those of the cCSl and eCS2 approximations. It is clear that for the wall- 
particle contact values the eCS3 approximation yields the best performance, 
while for the particle-particle contact values both the eCS2 and eCS3 are of 
comparable accuracy. A further feature to be pointed out is that the practical 
collapse on a common curve of the simulation data in Figs. 3 and 4 provide a 
posteriori support for the universality ansatz made in Eq. (16). 

As mentioned earlier, tlic;rc exist extra consistency conditions (sec for in- 
stance Ref. [12]) that one might use as well within our approach. Assuming 
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Z 

ij 

Fig. 3. Plot of the difference gij{ffij) — gfj^™"^{ffij) as a function of the parameter 

Zij = {(Jiaj / aij)M2/ Mz for hard spheres (d = 3) at a packing fraction jy = 0.49. 
The symbols arc simulation data for the single fluid (circle, Ref. [36]), three binary 
mixtures (squares, Rcf [49]) with 02/0-1 — 0.3 and Xi = 0.0625, 0.125, and 0.25, and 
a ternary mixture (triangles, Ref. [50]) with 02/01 = |, 03/01 = |, and xi = 0.1, 
X2 = 0.2. The lines are the PY approximation (- • • -), the SPT approximation 

(- • - •), the eCSl approximation (•••), the eCS2 approximation ( ), and the 

eCS3 approximation ( — ). 



that the ansatz (16) stiU holds, some of these conditions are related to the 
derivatives of Q with respect to z, namely 



dg{ri,z) 



dz 



3rj 



2(l-r?)2' 



d^g{v,z) 



dz^ 



3?7 



z=0 



d^g{rj,z) 



dz^ 



,PY 



0. 



(61) 



(62) 



(63) 



Interestingly enough, as shown by Eq. (52), condition (61) is already satisfied 
by our c3 approximation without having to be imposed. On the other hand, 
condition (63) implies ^3 = in the e3 scheme and thus it is only satisfied if 
Qs = gf^'^ , in which case we recover the SPT. Condition (62) is not fulfilled 
either by the SPT or by the e3 approximation (except for a particular ex- 
pression of Qs which is otherwise not very accurate). Thus, fulfilling the extra 
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Fig. 4. Plot of the difference g^j — g^f™"^ as a function of the parameter z-ujjll = 
(7jM2lM3 for hard spheres (d = 3) at a packing fraction r] — 0.4. The symbols are 
simulation data for a polydisperse mixture with a narrow top-hat distribution (open 
squares, Ref. [51]), a polydisperse mixture with a wide top-hat distribution (open 
circles, Ref. [51]), a polydisperse mixture with a Schulz distribution (open triangles, 
Ref. [51]), and a binary mixture (closed circles, Ref. [52]). The lines are the PY 
approximation (-•■-), the SPT approximation (- • - •), the cCSl approximation 
(•••), the eCS2 approximation ( ), and the eCS3 approximation ( — ). 

conditions (62) and (63) with a free requires either considering a higher 
order polynomial in z (in which case the consistency condition Z ~ can- 
not be satisfied for arbitrary mixtures, as discussed before) or not using the 
universality ansatz at all. In the first case, we have checked that a quartic or 
even a quintic polynomial docs not improve matters, whereas giving up the 
universality assumption increases significantly the number of parameters to 
be determined and seems not to be adequate in view of the behavior observed 
in the simulation data. 

An additional comment has to do with the restriction to d = 3 in this 
subsection. As noted before, the approximation el reduces to the exact result 
(3) for d = \. For d = 2, the approximation c2 already fulfills the condition 
Z = Z^ and so there is no real need to go further in that case. Since we 
have needed the approximation e3 to satisfy Z = for d = 3, it is tempting 
to speculate that a polynomial form for Q{z) of degree d could be found to 
be consistent with the condition Z = for d > 4. However, a detailed 
analysis shows that this is not the case for an arbitrary mixture, since the 
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number of conditions exceeds the number of unknowns, unless the universahty 
assumption is partially relaxed. 

As a final comment, let us stress that, although the discussion in this 
section has referred, for the sake of simplicity, to discrete mixtures, all the 
dependence on the details of the composition occurs through a finite number 
of moments, so that the results remain meaningful even for continuous poly- 
disperse mixtures [53]. In that case, instead of a set of mole fractions {xi} and 
a set of diameters {en}, one has to deal with a distribution function w{a) such 
that w{a)da is the fraction of particles with a diameter comprised between a 
and a + da. Therefore, the moments (1) are now defined as 



and with such a change the results we have derived for discrete mixtures also 
hold for polydisperse systems. 

2.3 Non- Additive Systems 

Non- additive hard-core mixtures, where the distance of closest approach be- 
tween particles of different species is no longer the arithmetic mean of the 
diameters of both particles, have received much less attention than additive 
mixtures, in spite of their in principle more versatility to deal with interesting 
aspects occurring in real systems (such as fluid-fluid phase separation) and 
of their potential use as reference systems in perturbation calculations on the 
thermodynamic and structural properties of, say, Lennard-Jones mixtures. 
Nevertheless, the study of non-additive systems goes back flfty years [54 -56] 
and is still a rapidly developing and challenging problem. 

As mentioned in the paper by Ballone et al. [57], where the relevant 
references may be found, experimental work on alloys, aqueous electrolyte 
solutions, and molten salts suggests that hetero-coordination and homo- 
coordination may be interpreted in terms of excluded volume effects due to 
non-additivity of the repulsive part of the intermolecular potential. In particu- 
lar, positive non-additivity leads naturally to demixing in HS mixtures, so that 
some of the experimental findings of phase separation in the above mentioned 
(real) systems may be accounted for by using a model of a binary mixture of 
(positive) non-additive HS. On the other hand, negative non-additivity seems 
to account well for chemical short-range order in amorphous and liquid binary 
mixtures with preferred hetero-coordination [58]. 

Some Preliminary Definitions 

Let us consider an A/'-component mixture of non-additive HS in d dimensions. 

In this case, aij = ■^{a.^ + aj)(l + Aij), where Aij > —1 is a symmetric 
matrix with zero diagonal elements {A^ = 0) that characterizes the degree 




(64) 
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of non-additivity of the interactions. If Aij > the non-additivity character 

of the ij interaction is said to be positive, while it is negative if Aij < 0. In 
the case of a binary mixture {N = 2), the only non-additivity parameter is 
A = Ai2 = A21. The virial EOS (2) remains being valid in the non-additive 
case. 

The contact values gij{(Jij) can be expanded in a power series in density 

as 

N N 

+ 0{p^). (65) 

k=l kj=l 

The coefficients Ck;ij , CkZ-ij , ■ • • are independent of the composition of the mix- 
ture, but they are in general complicated nonlinear functions of the diameters 
aij, (Jik, (Tjk, (Jke, ■ ■ ■ ■ Insertion of the expansion (65) into Eq. (2) yields the 
virial expansion of Z, namely 

00 

Z{p) = l + J2BniVdPr-' 

71=2 

JV JV 

= 'i- + VdP^ BijXiXj + {vdpf E BijkXiXjXk 

JV 

+{vdpf X] BijkeXiXjXkXe + 0{p'^). (66) 

i,j,k,e=l 

Note that, for further convenience, we have introduced the coefficients i3„ = 
v'^^"'~^^ Bn, where Bn are the usual virial coefficients [cf. Eq. (25)]. The 
composition-independent second, third, and fourth (barred) virial coefficients 
are given by 

Bij = 2''-'atj, (67) 
Bijk = {ck;ijcrfj + Cj-ikcyfk + Ci-jk<J%) , (68) 
2d-i 

Bijke = {cki;ij<yij + cji-ikafk + cu-jk<y% + CjkM<^ii + Cik,ji<y'je 

+Cir,ke4e) ■ (69) 



A Simple Proposal for the Equation of State of d-Dimensional 
Non- Additive Mixtures 

Our goal now is to generalize the el proposal given by Eq. (20) to the non- 
additive case [59]. We will not try to extend the e2 and e3 proposals, Eqs. (34) 
and (56), because of two reasons. First, given the inherent complexity of non- 
additive systems, we want to keep the approach as simple as possible. Second, 
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we are more interested in the EOS than in the contact values themselves and, 
as mentioned earlier, the cl proposal provides excellent EOS, at least in the 
additive case, despite the simplicity of the corresponding contact values. 

As the simplest possible extension, we impose again the point particle and 
equal size consistency conditions, Eqs. (12) and (13), and thus keep in this 
case also the ansatz (16) and the linear structure of Eq. (19). However, instead 
of using Eq. (15), we determine the parameters Zij as to reproduce Eq. (65) 
to first order in the density. The result is readily found to be [59] 



b2 J \ Md 



(70) 



Here 62 = 2'^"-'^ and 63 are the second and third virial coefficients for the 
single component fluid, as defined by Eq. (26). The proposal of Eq. (19) sup- 
plemented by Eq. (70) is, by construction, accurate for densities low enough 
as to justify the truncated approximation giji^ij) ~ 1 + VdPy^,^ XkCk-.g ■ On 
the other hand, the limitations of this truncated expansion for moderate and 
large densities may be compensated by the use of gs- When Eqs. (16), (19), 
and (70) are inserted into Eq. (2) one gets 



7( \ 1^ ^ b3MdB2--hB3 B3-MdB2 



Equation (71) is the sought generalization of Eq. (21) to non-additive hard- 
core systems. As in the additive case, the the density dependence in the EOS 
of the mixture is rather simple; Z{ri) — 1 is expressed as a linear combination 
of r]/{l — rf) and Zs{r}) — 1, with coefficients such that the second and third 
virial coefficients are reproduced. Again, Eq. (71) is bound to be accurate for 
sufficiently low densities, while the limitations of the truncated expansion for 
moderate and large densities are compensated by the use of the EOS of the 
pure fluid. 

The exact second virial coefficient B2 is known from Eq. (67). In principle, 
one should use the exact coefficients Ck-ij to compute B3. However, to the 
best of our knowledge they are only known for d < 3. Since our objective is to 
have a proposal which is explicit for any d, we can make use of a reasonable 
approximation for them [59], as described below. 



An Approximate Proposal for Ck- 



The values of the coefficients Ck;ij are exactly known for d = I and d = 3 and 
from these results one may approximate them in d dimensions as [59] 



'2 



Cfc;ij = (^k;ij + ( IT ~ ) -^r^(^i;jk(^j;ik^ (^2) 



where we have called 
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(^k;ij = (^ik + <^jk — (^ij (73) 

and it is understood that ak;ij > for all sets ijk. Clearly, ai-ij = ai. For a 
binary mixture Eq. (72) yields 

Cl;ll = {h/b2)(Tf, 

C2;ii = (2ai2 -aiY + (63/&2 - 1) (71 (2ai2 - ai)'^-\ (74) 

Cl;12 = af + (63/62 - 1) (2(712 - (7l)(7f/CTl2. 

Of course, Eqs. (72) and (74) reduce to the exact results for d = 1 (62 = 63 = 1) 
and for d = 3 (62 = 4, 63 = 10). 

The quantities ak;ij may be given a simple geometrical interpretation. 
Assume that we have three spheres of species and k aligned in the sequence 
ikj. In such a case, the distance of closest approach between the centers of 
spheres i and j is a^. + ajj^. If the sphere of species k were not there, that 
distance would of course be cr^ . Therefore a^-ij as given by Eq. (73) represents 
a kind of effective diameter of sphere k, as seen from the point of view of the 
interaction between spheres i and j. 

Inserting Eq. (72) into Eq. (70), one gets 

_fb3 A"VEfea^fc<^ , T,k^k(Tt7ijai-jk(Tj;ik 

~[h~V [ M, - ^ j + M^j • ^^^> 

It can be easily checked that in the additive case {crk;ij — * crfc), Eq. (75) reduces 
to Eq. (15). 

Equations (72) and (74) are restricted to the situation ak-ij > for any 
choice of i, j, and fc, i.e., 2cri2 > max(a"i, (T2) in the binary case. This excludes 
the possibility of dealing with mixtures with extremely high negative non- 
additivity in which one sphere of species k might "fit in" between two spheres 
of species i and j in contact. Since for d = 3 and N = 2 the coefficients Ck;ij 
are also known for such mixtures [60], we may extend our proposal to deal 
with these cases: 

Cl;ll = (63/62)0-^, 

C2;ll =^i + (63/62 - 1) 'yi^t\ (76) 
Cl;12 = (2(712 - + (63/62 - 1) CT2Crf /o-i2, 

where we have defined 

a2 = max (2iti2 - <7i , 0) . (77) 

With such an extension, we recover the exact values of Ck-ij for a binary 
mixture of hard spheres (d = 3), even if o"i > 2cri2 or (T2 > 2cri2. 

The EOS (71) becomes explicit when B3 is obtained from Eq. (68) by 
using the approximation (72). The resulting virial coefficient is the exact one 
for d = I and d ^ 3. For hard disks (d = 2), it turns out that the approximate 
third virial coefficient is practically indistinguishable from the exact one [59] . 
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Fig. 5. Plot of the compressibility factor versus the non-additivity parameter A for 

a symmetric binary mixture of non-additive hard spheres (d = 3) at rj — n/30 and 
two different compositions. The solid lines are our proposal, Eq. (71), with = Z^^, 
while the dashed lines arc Hamad's proposal (Refs. [61-63]). The symbols are results 
from Monte Carlo simulations (Refs. [64,65]). 

When the approximate B3 is used, Eq. (71) reduces to Eq. (21) in the additive 
case. 

From the comparison with simulation results, both for the compressibility 
factor and higher order virial coefficients, we find that the EOS (71) does a 
good job for non-additive mixtures, thus representing a reasonable compro- 
mise between simplicity and acciiracy, provided that Zg is accurate enough. 
This is ilhistrated in Fig. 5, where the proposal (71) with = Z^^ and a 
similar proposal by Hamad [61 63] are compared with simulation data [64,65] 
for some three-dimensional symmetric mixtures. A more extensive compari- 
son [59] shows that Eq. (71) seems to work better (especially as the density 
is increased) in the case of positive non-additivities, at least for d = 1, d = 2, 
and d = 3, but its performance is also reasonably good in highly asymmetric 
mixtures, even for negative A. Of course the full assessment of this proposal 
is still pending since it involves many facets (non-additivity parameters, size 
ratios, density, and composition). Without this full assessment and given its 
rather satisfactory performance so far, going beyond the approximation given 
by Eq. (19) (taking similar steps to the ones described in Subsections 2.1 and 
2.2 for additive systems) does not seem to be necessary at this stage, although 
it is in principle feasible. 
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2.4 Demixing 

Demixing is a common phase transition in fluid mixtures usually originated 
on the asymmetry of the interactions (e.g., their strength and/or range) be- 
tween the different components in the mixture. In the case of athcrmal systems 
such as HS mixtures in d dimensions, if fluid-fluid separation occurs, it would 
represent a neat example of an entropy-driven phase transition, i.e., a phase 
separation based only on the size asymmetry of the components. The exis- 
tence of demixing in binary additive three dimensional HS mixtures has been 
studied theoretically since decades, and the issue is still controversial. In this 
subsection wc will present our results following different but related routes 
that attempt to clarify some aspects of this problem. 

Binary Mixtures of Additive d-Dimensional Spheres (d = 3, d = 4 
and d = 5) 

Now wc look at the possible instability of a binary fluid mixture of HS of 
diameters cti and (72 (cri > (J2) in d dimensions by looking at the Helmholtz 
free energy per unit volume, /, which is given by 



where A, is the thermal de Broglie wavelength of species i. We locate the 
spinodals through the condition /11/22 — /12 — 0, with = d"^ f /dpidpj. Due 
to the spinodal instability, the mixture separates into two phases of different 
composition. The coexistence conditions are determined through the equality 
of the pressure p and the two chemical potentials ^1 and ^2 in both phases 
{jii = df /dpi), leading to binodal (or coexistence) curves. 

We begin with the case d = 3. It is well known that the BMCSL EOS, Eq. 
(11), does not lead to demixing. However, other EOS for HS mixtures have 
been shown to predict demixing [41,66], including the EOS that is obtained 
by truncating the virial series after a certain number of terms [67, 68] . In 
particular, it turns out that both Z = Z^^csi, Eq. (40), and Z = .^eCS2, 
Eq. (41). lead to demixing for certain values of the parameter 7 = a"2/o"i 
that measures the size asymmetry. The critical values of the pressure, the 
composition, and the packing fraction are presented in Table 2 for a few 
values of 7. 

As discussed earlier, the eCSl EOS and, to a lesser extent, the eCS2 EOS 
are both in reasonably good agreement with the available simulation results 
for the compressibility factor [18, 36, 45] and lead to the exact second and 
third virial coefficients but differ in the predictions for i?„ with n > 4. The 
scatter in the values for the critical constants shown in Table 2 is evident and 
so there is no indication as to whether one should prefer one equation over 
the other in connection with this problem. Notice, for instance, that the eCS2 




(78) 
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Table 2. Critical constants pcCri/ksT, xic, and rjc for different 7-values as obtained 
from the two extended CS ('(luatious (10) and (11). 





cCSl 


eCS2 


7 




T xic r/c 




T xic ric 


0.05 


3599 


0.0093 0.822 


1096 


0.0004 0.204 


0.1 


1307 


0.0203 0.757 


832.0 


0.0008 0.290 


0.2 


653.4 


0.0537 0.725 






0.3 


581.9 


0.0998 0.738 






0.4 


663.4 


0.1532 0.766 







does not predict demixing for 7 > 0.2, while both the values of the critical 
pressures and packing fractions for which it occurs according to the eCSl EOS 
suggest that the transition might be metastable with respect to a fluid-solid 
transition. 

Now we turn to the cases d = A and d = 5. Here we use the extended 
Luban-Michels equation (eLMl) described in Subsection 2.1 [see Eq. (21) 
and Table 1]. As seen in Fig. 6, the location of the critical point tends to go 
down and to the right in the 772 vs 771 plane as 7 decreases for d = 4 [69]. 
On the other hand, while it also tends to go down as 7 decreases if d = 5, 
its behavior in the r?2 vs r}i plane is rather more erratic in this case. Also, 
the value of the critical pressure pc (in units of ksT/af) is not a monotonic 
function of 7; its minimum value lies between 7=1/3 and 7 = 1/2 when 
d = 4, and it is around 7 = 3/5 for d = 5. This non- monotonic behavior is 
also observed for three-dimensional HS [66,68]. 

It is conceivable that the demixing transition in binary mixtures of hard 
hyperspheres in four and five dimensions described above may be metastable 
with respect to a fluid-solid transition, as it may also be the case of 3D HS. 
In fact, the value of the pressure at the freezing transition for the single 
component fluid is [31] pfa'^/ksT ~ 12.7 (d = 3), 11.5 (d = 4), and 12.2 
(d = 5), i.e., picr'^ /UbT does not change appreciably with the dimensionality 
but is clearly very small in comparison with the critical pressures pcfff/kBT 
we obtain for the mixture; for instance, pccrf/kBT ~ 600 (d = 3, 7 = 3/10), 
300 (d = 4, 7 = 1/3) and 123 (d = 5, 7 = 3/5). However, one should also 
bear in mind that, if the concentration xi of the bigger spheres decreases, the 
value of the pressure at which the solid-fluid transition in the mixture occurs 
in 3D is also considerably increased with respect to pf [cf. Fig. 6 of Ref. [66]]. 
Thus, for concentrations Xi ~ 0.01 corresponding to the critical point of the 
fluid-fluid transition, the maximum pressure of the fluid phase greatly exceeds 
Pf. If a similar trend with composition also holds in 4D and 5D, and given that 
the critical pressures become smaller as the dimensionality d is increased, it is 
not clear whether the competition between the fluid-solid and the fluid-fluid 
transitions in these dimensionalities will always be won by the former. The 
point clearly deserves further investigation. 
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Fig. 6. Spinodal curves (upper panels: lines) and binodal curves (upper panels: open 

symbols; lower panels: linos) in a 4D system (left panels) and in a 5D system (right 
panels). The closed symbols are the critical consolute points. 



An interesting feature must be mentioned. There is a remarkable similarity 

between the binodal curves represented in the paf rji and in the /i; rji planes 
[69]. By eliminating rji as if it were a parameter, one can represent the binodal 
curves in a /Zi vs paf plane. Provided the origin of the chemical potentials 
is such as to make = (j,;, the binodals in the jii-paf plane practically 
collapse into a single curve (which is in fact almost a straight line) for each 
dimensionality (d = 3, d = 4, and = 5) [69]. A closer analysis of this 
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phenomenon shows, however, that it is mainly due to the influence on /ij 
of terms which arc quantitatively dominant but otherwise irrelevant to the 
coexistence conditions. 

Binary Mixtures of Non- Additive Hard Hyperspheres in the Limit 
of High Dimensionality 

Let us now consider a binary mixture of non-additive HS of diameters ai 
and (72 in d dimensions. Thus in this case ai2 = ^{cti + 0-2) (1 + A) where as 
before A may be either positive or negative. Further assume (something that 
will become exact in the limit d ^ 00 [70]) that the EOS of the mixture is 
described by the second virial coefficient only, namely 

p = pkBT[l + B2{xi)p\, (79) 

where, according to Eq. (67), 

B^ixi) = Vd2'^-^ {xiaf + xl4 + 2x1x2(7(2) . (80) 

The Helmholtz free energy per unit volume is given by f/phsT = — 1 + 
S^i=i 111 {xipXf) +B2P, where Eq. (78) has been used. The Gibbs free energy 
per particle is 

2 

9 = U+p)Ip = X^a;i In {xipXf) + 2B2{xi)p, (81) 

where without loss of generality we have set ksT = 1. Given a size ratio 7, a 
value of Z\, and a dimensionality d, the consolutc critical point {xicPc) is the 
solution to (d'^g/dxf)^ = (d'^g/dxf)^ = 0, provided of course it exists. Then, 
one can get the critical density Pc from Eq. (79). 
We now introduce the scaled quantities [71] 

p = 2^-'^Vdd-'^paf/kBT, u = rf-^Bap. (82) 

Consequently, Eqs. (79) and (81) can be rewritten as 

p = u{u + d-^) /B2, (83) 

2 

g = '^Xiln (xiAi) + In (yl^u/i^z) + 2du, (84) 

i=l 

where B2 = B2/2''-'^Vd(7f, A,, = (Xi/aiY, and Ad = d/2'^-'^Vd. Next we take 
the limit d ^ 00 and assume that the volume ratio 7 = 7'' is kept fixed and 
that there is a (slight) non-additivity A = d~^A such that the scaled non- 
additivity parameter A is also kept fixed in this limit. Thus, the second virial 
coefficient can be approximated by 
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(85) 

with 

J = i (In 7)' + 2 A (86) 

Let us remark that, in order to find a consolutc critical point, it is essential 
to keep the term of order ii A < 0. The EOS (83) can then be inverted 
to yield 



2 V ^2 / 

(87) 

In turn, the Gibbs free energy (84) becomes 
5 = 5(")rf + 5(i)+0(d-i), 

while the chemical potentials /xi = g+X2 {dg/dxij^ and /i2 = g—x\ {dg/dxi)^ 
are given by 

M« = In (a,x,A,^I^) - 1/^ + {x2/x,){W"b'^'Ib'^\ ^"'^ 



where (12 is obtained from /^i by the changes a;i X2, vli — > 712/7, 7 — > I/7, 

p ^P7, B2 ^ ^2/7- 

The coordinates of the critical point are readily found to be 

73/4 (1 + 7^/^)^ 

^^^ = TTW^' 47P • 

Note that xi,; is independent of Z\. The coexistence curve, which has to be 
obtained numerically, follows from the conditions iJii\xA,p) — \^i\xB,p) 
(i = 1, 2) where xi = xa and xi = xb are the mole fractions of the co- 
existing phases. Once the critical consolute point has been identified in the 
pressure/concentration plane, we can obtain the critical density. The domi- 
nant behaviors of B2 and u at the critical point are 

' ™ (l_;^i/4 + ;^i/2)2' 2(1-71/4 + 71/2) J- ^^'^ 

Hence, the critical density readily follows after substitution in the scaling 
relation given in Eq. (82). It is also convenient to consider the scaled version 
T] = d~^2'^ri of the packing fraction 77 = vapaf {xi + X2'y)- At the critical point, 
it takes the nice expression 




^1 



Fig. 7. Binodal curves in the planes p vs x\ and 77 vs x\ corresponding to 7 = 0.01 
and A = -0.1, Z\ = 0, and Zi = 0.1. 



Vc = J ^. (92) 

The previous results clearly indicate that a demixing transition is possible 
not only for additive or positively non-additive mixtures but even for negative 
non-additivities. The only requirement is J > 0, i.e., A > — | (In 7)^ or, 
equivalently, A > — i (In 7)^. Figure 7 shows the binodal curves corresponding 
to 7 = 0.01 and zi = — 0.1 (negative non-additivity), A = (additivity), and 
Z\ = 0.1 (positive non-additivity). 

While the high dimensionality limit has allowed us to address the prob- 
lem in a mathematically simple and clear-cut way, the possibility of demixing 
with negative non-additivity is not an artifact of that limit. As said before, 
demixing is known to occur for positive non-additive binary mixtures of HS 
in three dimensions and there is compelling evidence on the existence of this 
phenomenon in the additive case, at least in the metastable fluid region. Even 
though in a three-dimensional mixture the EOS is certainly more complicated 
than Eq. (79) and the demixing transition that we have just discussed for neg- 
ative non-additivity is possibly metastable with respect to the freezing transi- 
tion, the main effects at work (namely the competition between depletion due 
to size asymmetry and hetero-coordination due to negative non-additivity) 
are also present. In fact, it is interesting to point out that Roth et al. [72], 
using the approximation of an effective single component fluid with pair inter- 
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actions to describe a binary mixture of non-additive 3D HS and employing an 
empirical rule based on the effective second virial coefficient, have also sug- 
gested that demixing is possible for small negative non-additivity and high 
size asymmetry. Our exact results lend support to this suggestion and con- 
firm that, in some cases, the limit d — > oo highlights features already present 
in real systems. 



3 The Rational Function Approximation (RFA) Method 
for the Structure of Hard-Sphere Fluids 

The RDF g{r) and its close relative the (static) structure factor S{q) are the 
basic quantities used to discuss the structure of a single component fluid [f-4]. 
The latter quantity is defined as 

S{q) = l+ph{q), (93) 

where 

h{q) = I dre-''i '"/i(r) (94) 



is the Fourier transform of the total correlation function h{r) = g{r) — l,\ being 
the imaginary unit. An important related quantity is the direct correlation 
function c(r) , which is defined in Fourier space through the Ornstein-Zernike 
(OZ) relation [1-4] 

1 -I- ph{q) 

where c(g) is the Fourier transform of c(r) 

The usual approach to obtain g{r) is through one of the integral equa- 
tion theories, where the OZ equation is complemented by a closure relation 
between c(r) and h{r) [1]. However, apart from requiring in general hard nu- 
merical labor, a disappointing aspect is that the substitution of the (necessar- 
ily) approximate values of g{r) obtained from them in the (exact) statistical 
mechanical formulae may lead to the thermodynamic inconsistency problem. 

The two basic routes to obtain the EOS of a single component fluid of HS 
are the virial route, Eq. (14), and the compressibility route 

Xs = ksT (^) = [1 - pc{0)]-' = S{0) 



dp 



T 



1 + drr''-^h{r). (96) 

Jo 



Thermodynamic consistency implies that 



X;HV) = ^ivZsiv)], (97) 



Alternative Approaches to Hard-Sphere Liquids 



31 



but, in general, this condition is not satisfied by an approximate RDF. In 
the case of a HS mixture, the virial route is given by Eq. (2), while the 
compressibility route is indicated below [cf. Eq. (145)]. 

In this section we describe the RFA method, which is an alternative to 
the integral equation approach and in particular leads by construction to 
thermodynamic consistency. 



3.1 The Single Component HS Fluid 

We begin with the case of a single component fluid of HS of diameter u. 
The following presentation is equivalent to the one given in Rcfs. [73,74], 
where all details can be found, but more suitable than the former for direct 
generalization to the case of mixtures. 

The starting point will be the Laplace transform 

POO 

G{s) = / dre-*V5(r) (98) 
Jq 

and the auxiliary function ^{s) defined through 

G(s) = ^[p + e-!Z'(5)]-\ (99) 

The choice of G{s) as the Laplace transform of rg{r) and the definition of 
!Z'(s) from Eq. (99) are suggested by the exact form of g{r) to first order in 
density [73]. 

Since g{r) = for r < tr while g{<J^) = finite, one has 

g{r) = 0(r - a) [g{a+) + g'{<y+){r - a) + • • • ] , (100) 

where g'{r) = dg{r)/dr. This property imposes a constraint on the large s 
behavior of G(s), namely 

e-^'sGis) = ag{a+) + [gia+) + ag'{a+)] s-' + ©(s'^). (101) 

Therefore, lims^ao e^" sG{s) = ag{a^) = finite or, equivalently, 

lim s-'^^(s) = = finite. (102) 

On the other hand, according to Eq. (96) with d = 3, 

d (°° 

Xs = 1 - 24r7cr"^ lim — / dr e"^''r \g{r) - 1] 
ds Jo 

= 1 - 24r/a-^ lini ^ [G(s) - s'^] . (103) 

Since the (reduced) isothermal compressibility Xs is also finite, one has 

dr [5(^) ^ 1] = finite, so that the weaker condition drr \g(r) — 1] = 
lims^o[G(s) — s~^] = finite must hold. This in turn implies 
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^{s) = -p + pas - ^pah^ + Qpa^ + i-^ - (J^pa^ + ^) a/ + O(s^). 

(104) 

First-Order Approximation (PY Solution) 

An interesting aspect to be remarked is that the minimal input we have just 
described on the physical requirements related to the structure and thermo- 
dynamics of the system is enough to determine the small and large s limits of 
!Z'(s), Eqs. (102) and (104), respectively. While infinite choices for '?'(s) would 
comply with such limits, a particularly simple form is a rational function. In 
particular, the rational function having the least number of coefBcients to be 
determined is 

'^(^^ = mnWs ' ^'''^ 

where one of the coefficients can be given an arbitrary non-zero value. We 
choose E'-^^ = 1. With such a choice and in view of Eq. (104), one finds 
£(0) = -pL(0), £(1) = -p(L(i) - aiW), E(^) = p{aL^^) - ia^iW), and 

^^°^ = 2-(^> (106) 

Upon substitution of these results into Eqs. (99) and (105), we get 

^^^^ " 27rs2 1 - p [^2{(Ts)a^m + Vi{<^s)a^LW] ' ^^^^^ 



where 



In particular. 



\m=0 



^^{x) ^ ( ^ - e-- ) . (109) 



1-e-^ _ l-x-c""' _ 1 - .)• + .i:-'/2 - o 



— X 



^o{x) = , tpiix) = 2 ' ¥'2(a;) 

(110) 

Note that lim^^o ^n{x) = {-!)'"' /{n + 1)!. 

It is remarkable that Eq. (108), which has been derived here as the sim- 
plest rational form for iP^{s) complying with the requirements (102) and (104), 
coincides with the solution to the PY closure, c(r) = for r > cr, of the OZ 
equation [29]. Application of Eq. (102) yields the PY contact value gf^ and 
compressibility factor Zf^ shown in Table 1. Analogously, Eq. (103) yields 



(1 + 277)2 



xr^ = irriz^- (111) 



It can be easily checked that the thermodynamic relation (97) is not satisfied 
by the PY theory. 
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In the spirit of the RFA, the simplest extension of the rational approximation 
(105) involves two new terms, namely as* in the numerator and i^^^s^ in the 
denominator, both of them necessary in order to satisfy Eq. (102). Such an 
addition leads to 



,£;(3)53 



i(0) +i(l)5 + i(2)32 



Applying Eq. (104), it is possible to express 
L(^) in terms of a and L^^). This leads to 



(112) 



E(^\ L(0), and 



G{s) 



where 



(113) 



= 277 



L(^) = 27rc7- 



1 + 2t/ 



127] 
1 — T/ 



+ 





a 


L(2)\ 


1 - 


7] a 


" CT2 j' 


1 + 
1 - 


— a 
V 


L(2) 

-37? 



(114) 



(115) 



(l-r?)2 ■ 1-n 

Thus far, irrespective of the values of the coefficients ^^2) and a, the condi- 
tions lims^oc c^'^'SG'(.s) = finite and lims^o[G(s) — s~2j — finite arc satisfied. 
Of course, if L'^2) = q, = q, one recovers the PY approximation. More gen- 
erally, we may determine these coefficients by prescribing the compressibility 
factor Zs (or equivalently the contact value gs) and then, in order to ensure 
thermodynamic consistency, compute from it the isothermal compressibility 
Xs by means of Eq. (97). From Eqs. (102) and (103) one gets 



(l(0)) 



L^^^ = 2Traag, 
12r) a 



iZTi a „a\ izr? ah 

1-- — l-h2- -h — - — : 

1 — ryaV cr/ tt o" 



1277 aL(2) 



(116) 
(117) 



Clearly, upon substitution of Eqs. (114) and (116) into Eq. (117) a quadratic 
algebraic equation for a is obtained. The physical root is 



a = — 



12r/(l + 2T])Ei 



(l-7?)2+367?[l + r?-Z,(l-7?)]£;4' 



where 



367? {Z, - \) 



z. - zr 



,/PY 



1/2 ~ 



(118) 



(119) 



The other root must bo discarded because it corresponds to a negative value 
of a, which, according to Eq. (116), yields a negative value of L(2)_ xhis would 
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imply the existence of a positive real value of s at which G{s) = [73,74], 
which is not compatible with a positive definite RDF. However, according to 
the form of Eq. (119) it may well happen that, once Zs has been chosen, there 
exists a certain packing fraction r]g above which a is no longer positive. This 
may be interpreted as an indication that, at the packing fraction rjg whore a 
vanishes, the system ceases to be a fluid and a glass transition in the HS fluid 
occurs [74-76]. 

Expanding (113) in powers of s and using Eq. (101) one can obtain the 
derivatives of the RDF at r = cr+ [77]. In particular, the flrst derivative is 



(120) 



which may have some use in connection with perturbation theory [15]. 

It is worthwhile to point out that the structure implied by Eq. (113) coin- 
cides in this single component case with the solution of the Generalized Mean 
Spherical Approximation (GMSA) [78], where the OZ relation is solved under 
the ansatz that the direct correlation function has a Yukawa form outside the 
core. 

For a given Zg, once G(s) has been determined, inverse Laplace trans- 
formation yields rg{r). First, note that Eq. (99) can be formally rewritten 
as 

oo 

^(^) = E/^""' [-•^^(sr^e-"^^ (121) 

n— 1 

Thus, the RDF is then given by 



-H oo 

with {x) denoting the Heaviside step function and 

^„(r) = -£-i{s[-!Z^(s)]-"}, (123) 

denoting the inverse Laplace transform. Explicitly, using the residue the- 
orem, 

4 n (i) 

V'nW = -Ee'"'E7 ITT nr^""™' (124) 

, ^An — myAm—ly. 

1=1 TO=1 ^ ' ^ ' 



where 



1 \ m— 1 

d 



aW = lim ( - ) s[-^ {s) lis - s,)]"" , (125) 



Si (i = 1, ... ,4) being the poles of l/!f'(s), i.e., the roots of E^^'^ + E^'^h + 
E^'^\s'^ + ij('^)s'^ + as** = 0. Explicit expressions of g{r) up to the second 
coordination shell a < r < 3a can be found in Ref. [79]. 
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On the other hand, the static structure factor S{q) [ci. Eq. (93)] and the 
Fourier transform h{q) may be related to G{s) by noting that 



h{q) 



Aw 
Q 



f 

Jo 



drrsm{qr)h{r) = — 27r 



Gis)-G{-s) 



(126) 



Therefore, the basic structural quantities of the single component HS fluid, 
namely the RDF and the static structure factor, may be analytically deter- 
mined within the RFA method once the compressibility factor Zg, or equiva- 
lently the contact value g^, is specified. In Fig. 8 we compare simulation data 
of g{r) for a density pa^ = 0.9 [80] with the RFA prediction and a recent 
approach by Trokhymchuk et al. [81], where Zs = Z^^ [cf. Table 1] and the 
associated compressibility 



X: 



cs 



1 + 477 + 4?7^ - 4773 



(127) 



are taken in both cases. Both theories are rather accurate, but the RFA cap- 
tures better the maxima and minima of g{r) [82]. 

It is also possible to obtain within the RFA method the direct correlation 
function c(r). Using Eqs. (95) and (126), and applying the residue theorem, 
one gets, after some algebra, 



5^ 

"Ho 




r 

Fig. 8. Radial distribution function of a single component HS fluid for pa^ = 0.9. 
The solid lines represent simulation data [80]. The dashed lines represent the results 

of the approach of Rcf. [81], while the dotted lines refer to those of the RFA method. 
The inset shows the oscillations of g{r) in more detail. 
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c(r) = [k+ + K_ + +Ko + Kir + K^r^ 0(1 -r)+K- 



where 



r r r j r 

(128) 

K = ^'^12ar?L(2)/7r+ 1 - 12a(l + 2a)r?/(l - r?), (129) 



^± = 4^^5(r-^{^ + + ^"^'^l ± [2 + r? + 2a(l + 2r?)] k 

+ (1 - ry) [k^ - ?7 (12 + (k ± 6)«;)] L*^)/,^! |l2r7 [1 + 2(1 + 3a)r]] 
±6r] [3r] - 2a{l - 4?])] k - 6r]{l + 2a)(l - t])k^ - (1 - rifn^ian T 1) 
+6r?(l - 77) [k^ - (12 + (k ± 6)k)] L^^^ /7r| , (130) 

i^-i = - + i^+e'= + if_e-« + ifo + i^i + , (131) 



1 + 2 (1 + 3a) r? - 67? (1 - r?) L^^)/^ ' ^ 
aK (1 — vif 



(132) 



X (7 + r? + 6a (2 + r?))] + 12^ (2 + r/) (1 - r?)2L(2)'/7r2} ,(133) 

^3 = \k^, (134) 
/<:=-(/<:+ + i^_+i^_i). (135) 

In Eqs. (129)-(135) we have taken a = 1 as the length unit. Note that Eq. 
(135) guarantees that c(0) = finite, while Eq. (131) yields c(cr+) — c(cr~) = 
L^'^^ /2'Ka = g{<y^)- The latter equation proves the continuity of the indirect 
correlation function 7(r) = h{r) — c(r) at r = a. With the above results, 
Eqs. (122) and (128), one may immediately write the function 7(r). Finally, 
we note that the bridge function B{r) is linked to 7(r) and to the cavity 
(or backgroimd) function y{r) = e'^^'^^^'^^'^ g{r), where (f>{r) is the interaction 
potential, through 

B(r)=lny(r)-7(r), (136) 

and so, within the RFA method, the bridge function is also completely speci- 
fied analytic;ally for r > a once Zg is prescribed. 

If one wants to have B{r) also for < r < cr, then an expression for the 
cavity function is required in that region. Here we propose such an expression 
using a limited number of constraints. First, since the cavity function and its 
first derivative are continuous at r = a, we have 
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0.0 0.2 0.4 0.6 0.8 1.0 

r 

Fig. 9. Cavity function of a single component HS fluid in the overlap region for 
pa^ = 0.3, 0.5, and 0.7. The sohd lines represent our proposal (140) with Zg = Z^^, 
while the symbols represent Monte Carlo simulation results [84]. 



v'(l) 1 

where Eqs. (116) and (120) have been used and again u = 1 has been taken. 
Next, we consider the following exact zero-separation theorems [83]: 

z (n') - 1 

In 2/(0) = ZM - 1 + / dr,'^^, , (138) 

Jo V 

^ = -^vy{i). (139) 

The four conditions (137)-(139) can be enforced by assuming a cubic poly- 
nomial form for In y{r) inside the core, namely 

y{r) = exp (Yq + Y^r + Y2r^ + Ysr^) , (0 < r < 1), (140) 

where 

Yo = ZM - 1 + r (141) 
Jo V 

Y, = -6^(1), (142) 
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Fig. 10. Parametric plot of the bridge function B{r) versus the indirect correlation 

function 7(r). The dashed line refers to the RFA for r) — 0.3, while the solid line 
refers to the RFA for rj — 0.49. In each case, the branch of the curve to the right 
of the circle corresponds to r < 1, while that to the left corresponds to r > 1. For 
comparison, the PY closure B{r) = ln[l + 7(r)] — 7(r) is also plotted (dash-dotted 
line) . 



y2 = 31ny(l)-ffl-3yo-2Fi, (143) 
Ys = -2lny{l) + ^ + 2Yo + Y,. (144) 

y(i) 

The proposal (140) is compared with available Monte Carlo data [84] in Fig. 
9, where an excellent agreement can be; observed. 

Once the cavity function y{r) provided by the RFA method is comple- 
mented by (140), the bridge function B{r) can be obtained at any distance. 
Figure 10 presents a parametric plot of the bridge function versus the indirect 
correlation function as given by the RFA method for two different packing 
fractions, as well as the result associated with the PY closure. The fact that 
one gets a smooth curve means that within the RFA the oscillations in j{r) 
are highly correlated to those of B{r). Further, the effective closure relation in 
the RFA turns out to be density dependent, in contrast with what occurs for 
the PY theory. Note that the absolute value |-B(r)| for a given value of 7(r) 
is smaller in the RFA than the PY value and that the RFA and PY curves 
become paradoxically closer for larger densities. Since the PY theory is known 
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to yield rather poor values of the cavity function inside the core [85,86], it 
seems likely that the present differences may represent yet another manifes- 
tation of the superiority of the RFA method, a point that certainly deserves 
to be further explored. 

3.2 The Multicomponent HS Fluid 

The method outlined in the preceding subsection will be now extended to an 
A''-component mixture of additive HS. Note that in a multicomponent system 
the isothermal compressibility x is given by 

N 

= 1 - p ^ XiXjCij{0), (145) 

where Cij{q) is the Fourier transform of the direct correlation function Cij{r), 
which is defined by the OZ equation 

N 

hij{q) = Cij{q) + ^ Pkhik{q)ckj{q), (146) 
fe=i 

where hij{r) = gij{r) — 1. Equations (145) and (146) are the multicompo- 
nent extensions of Eqs. (96) and (95), respectively. Introducing the quantities 
hij{q) = y/piPjhij{q) and Cij{q) = y/PiPj Cij{q), the OZ relation (146) be- 
comes, in matrix notation, 

c{q) = h{q)-[\ + h{q)]-\ (147) 

where I is the N x N identity matrix. Thus, Eq. (145) can be rewritten as 

JV JV 

Similarly to what we did in the single component case, we introduce the 
Laplace transforms of rgij{r): 

/.oo 

Gij{s) = / dre— r5i,(r). (149) 
Jo 

The counterparts of Eqs. (100) and (101) are 



l + h(0) 



(148) 
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e--«sGi,-(5) = a^^QiMj) + [9iMt^ + <^iAMj)\ + ^(«"')- (151) 

Moreover, the condition of a finite compressibility implies tliat hij{0) — finite. 
As a consequence, for small s, 

s'G.^is) = l + HS^s' + HiPs^ + ... (152) 
with H^j^ = finite and H-^^ = —hij{0)/ATr = finite, where 

dri-rrrhijir). (153) 



1 /•<- 
n\ Jo 



Wc arc now in the position to generalize the approximation (113) to the 
A^-component case [87]. While such a generalization may be approached in a 
variety of ways, two motivations are apparent. On the one hand, we want to 
recover the PY result [8] as a particular case in much the same fashion as in 
the single component system. On the other hand, we want to maintain the 
development as simple as possible. Taking all of this into account, we propose 

where L(s) and A(s) are the matrices 

L,,is)=4f+Ll]^s + 4fs\ (155) 

Aij{s) = Pi ip2{(Tis)afL\f + ipi{ais)a^L\f + ipo{ais)aiL\f , (156) 

the functions (/3„(.t) being defined by Eq. (109). We note that, by construc- 
tion, Eq. (154) complies with the requirement lims_»oo e'^'^^sG^ (s) = finite. 
Further, in view of Eq. (152), the coefficients of s° and s in the power series 
expansion of s^Gij{s) must be 1 and 0, respectively. This yields 2N'^ condi- 
tions that allow us to express 1^°^ and L^^) in terms of L^^^ and a. The solution 
is [87] 

N 

Lf) = ^1 + §2<yj + 2^2a - ^ PkOkL^k] ' (157) 

fc=i 



r(l) 



1 1 ^ 

= -dicJij + ■^^2cri(7j + (o^i + '&2(Ti)a - -'didi'^pkCTkL^^] , (158) 



fe=i 



where -di = 2^/(1 - rj) and ■i?2 = 6^(1^2 /Ma)?] /{I - -qf. 

In parallel with the development of the single component case, L*^^) and a 

(2) 

can be chosen arbitrarily. Again, the choice L-^- = a = gives the PY solution 
[8,88]. Since we want to go beyond this approximation, we will determine 
those coefficients by taking prescribed values for gjj((Ty ), which in turn, via 
Eq. (2), give the EOS of the mixture. This also leads to the required value of 
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X ^ = d{pZ)/dp, thus making the theory thermodynamically consistent. In 
particular, according to Eq. (151), 

4f =2naaijgij{a±). (159) 

The condition related to x is more involved. Making use of Eq. (152), one can 
get hij{0) = —AttHIP in terms of L^^^ and a and then insert it into Eq. (148). 

(2) 

Finally, elimination of Llj in favor of a from Eq. (159) produces an algebraic 
equation of degree 2A'", whose physical root is determined by the requirement 
that Gij{s) is positive definite for positive real s. It turns out that the physical 
solution corresponds to the smallest of the real roots. Once a is known, upon 
substitution into Eqs. (154), (157), (158), and (159), the scheme is complete. 
Also, using Eq. (151), one can easily derive the result 



1 


r(l) r(2) 






2'Kaaij 









(160) 



It is straightforward to check that the results of the preceding subsection are 
recovered by setting ai = cr, regardless of the values of the mole factions. 

Once Gij{s) has been determined, inverse Laplace transformation directly 
yields rgij{r). Although in principle this can be done analytically, it is more 
practical to use one of the efficient methods discussed by Abate and Whitt [89] 
to numerically invert Laplace transforms [90]. 

In Fig. 11 we present a comparison between the results of the RFA method 
with the FY theory and simulation data [50] for the RDF of a ternary mixture. 
In the case of the RFA, we have used the eCS2 contact values and the cor- 
responding isothermal compressibility. The improvement of the RFA over the 
FY prediction, particularly in the region near contact, is noticeable. Although 
the RFA accounts nicely for the observed oscillations, it seems to somewhat 
overestimate the depth of the first minimum. 

Explicit knowledge of Gij (s) also allows us to determine the Fourier trans- 
form hij {q) through the relation 



(161) 

if 

The structure factor Sij{q) may be expressed in terms of hij{q) as [4] 

Sij (q) = x^Sij + px^Xjhij (q). (162) 

In the particular case of a binary mixture, rather than the individual structure 
factors Sij{q), it is some combination of them which may be easily associated 
with fluctuations of the thermodynamic variables [91,92]. Specifically, the 
quantities [4] 

Snn{q) = Sn{q) + ^22(9) + 2Si2{q), (163) 
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Fig. 11. Radial distribution functions gij{r) for a ternary mixture with diameters 
CTi = 1, (72 = 2, and (73 = 3 at a packing fraction 77 = 0.49 with mole fractions 
XI = 0.7, X2 = 0.2, and xa = 0.1. The circles are simulation results [50], the solid 
lines are the RFA predictions, and the dotted lines are the PY predictions. 

Snc{q) = X2Sii{q) - XiS22{q) + " Xi)Si2{q), (164) 

S,,{q) = xlS^M + xlS22{q) - 2x,X2SM (165) 

are sometimes required. ^ 

After replacement of hij{q) = ^/PiPjhij{q) in Eq. (147), one easily gets 

Cij{q). Subsequent inverse Fourier transformation yields Cij{r). The result 
gives Cij (r) for r > cTy as the superposition of N Yukawas [93] , namely 
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N 



c,^{r)=Y,K\f-^, (166) 



-K£r 



where q = ±iK£ with £ = 1, ... ,7V are the zeros of det 



and the 



amplitudes K^' are obtained by applying the residue theorem as 

^If = ^ lim c„(g)(g-i«,). (167) 

The indirect correlation functions 7^ (r) = hij (r) — Cij (r) readily follow 
from the previous results for the RDF and direct correlation functions. Finally, 
in this case the bridge functions Bij{r) for r > Uij are linked to gij{r) and 
Cij{r) through 

i?,,(r)=ln.g,,(r)-7,,(r) (168) 

and so once more we have a full set of analytical results for the structural 
properties of a multicomponent fluid mixture of HS once the contact values 
9ij{^ij) specified. 



4 Other Related Systems 

The philosophy behind the RFA method to derive the structural properties 

of three-dimensional HS systems can be adapted to deal with other related 
systems. The main common features of the RFA can be summarized as follows. 
First, one chooses to represent the RDF in Laplace space. Next, using as a 
guide the low-density form of the Laplace transform, an auxiliary function is 
defined which is approximated by a rational or a rational-like form. Finally, 
the coefficients are determined by imposing some basic consistency conditions. 
In this section we consider the cases of sticky-hard-sphere, square-well, and 
hard-disk fluids. In the two former cases the RFA program is followed quite 
literally, while in the latter case it is done more indirectly through the RFA 
method as applied to hard rods {d = 1) and hard spheres (rf = 3). 

4.1 Sticky Hard Spheres 

The sticky-hard-sphere (SHS) fluid model has received a lot of attention since 
it was first introduced by Baxter in 1968 [94] and later extended to multi- 
component mixtures by Perram and Smith [95] and, independently, by Bar- 
boy [96] . In this model, the molecular interaction may be defined via square- 
well (SW) potentials of infinite depth and vanishing width, thus embodying 
the two essential characteristics of real molecular interactions, namely a harsh 
repulsion and an attractive part. In spite of their known shortcomings [97], an 
important feature of SHS systems is that they allow for an exact solution of the 
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OZ equation in the PY approximation [94,95]. Furthermore, they are thought 

to be appropriate for describing structural properties of colloidal systems, mi- 
celles, and microemulsions, as well as some aspects of gas-liquid equilibrium, 
ionic fluids and mixtures, solvent mediated forces, adsorption phenomena, 
polydispcrsc systems, and fluids containing chainlikc molecules [98-102]. 

Let us consider an A^-component mixture of spherical particles interacting 
according to the SW potential 



00, 



(f>ij{r) = 



r < aij, 

<r < Rij, 



0, r>R, 



(169) 



As in the case of additive HS, aij = {ai + aj)/2 is the distance between the 
centers of a sphere of species i and a sphere of species j at contact. In addition, 
Cij is the well depth and Rtj — a^j indicates the well width. We now take the 
SHS limit [94], namely 



Ri 



oo, 



12 Ri. 



finite. 



(170) 



where the are monotonically increasing functions of the temperature T and 
their inverses measure the degree of "adhesiveness" of the interacting spheres i 
and j. Even without strictly taking the mathematical limits (170), short-range 
SW fluids can be well described in practice by the SHS model [103]. 
The virial EOS for the SHS mixture is given by 



1 " /■ d 

i.7=l 



-<l>i,{r)lkBT 



N 



12t.. 



(171) 



where yij{r) = gij{r)e'^*^^^^/^'^'^ is the cavity function and ylj{r) = dyij{r)/dr. 
Since yij{r) must be continuous, it follows that 



9ij{r) =yij{r) 



(172) 



The case of a HS system is recovered by taking the limit of vanishing adhe- 
siveness Tj"^ ^ 0, in which case Eq. (171) reduces to the three-dimensional 
version of Eq. (2). On the other hand, the compressibility EOS, Eq. (145), is 
valid for any interaction potential, including SHS. 

As in the case of HS, it is convenient to define the Laplace transform (149). 
The condition yij{aij) = finite translates into the following large s behavior 
of Gij{s): 



e^^^^Gijis) 



( 1 



+ 0(s-'), 



(173) 
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which differs from (151): while e<"«*Gij(s) ~ s'^ for HS, e'"«*Gij(s) ~ 5° for 
SHS. However, the small s behavior is still given by Eq. (152), as a consequence 
of the condition = finite. 

The RFA proposal for SHS mixtures [104] keeps the form (154), except 
that now 

L.,is)^Lf;^+Lfs + Lfs^ + L^s^ (174) 

Aij{s) = Pi [p2{<^is)afLf + <pi(ais)a2Lg) + ^o{cTis)aiLlf - e-'^X^f 

(175) 

instead of Eqs. (155) and (156). By construction, Eqs. (154), (174), and (175) 
comply with the requirement lims^oo e'^*^^G'ij(s) = finite. Further, in view 
of Eq. (152), the coefficients of s° and s in the power series expansion of 
s'^Gij{s) must be 1 and 0, respectively. This yields 2N'^ conditions that allow 
us to express L^^^ and L^^^ in terms of L^^^, L^^^, and a as [104] 

N N 

fe=i fe=i 



1 1 ^ 

fc=i 

1 ^ 

-^{^i+'d20i)Y.PkOkLtl (177) 



where 'di and ^2 are defined below Eq. (158). Wc have the freedom to choose 
L^'^) and a, but L^^^ is constrained by the condition (173), i.e., the ratio 
between the first and second terms in the expansion of e^'i^Gij{s) for large s 
must be exactly equal to aij/12Tij. 

First-Order Approximation (PY Solution) 

The simplest approximation consists of making o: = 0. In view of the condition 
e'^*'^Gij{s) ~ s° for large s, this implies ig' = 0. In that case, the large s 
behavior that follows from Eq. (154) is 



27re'^«^Gy(s)=Lg) 

where 



L<]> + (L<=.^D)^^ 



-i+0(s-2), (178) 



D.,^P.{^a^L'ff-a.L^+Lfj. (179) 
Comparison with Eq. (173) yields 



ViMj) = ^Lf^, (180) 
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1 9_. . r (2) N 

^^=^+i:^D,,. (181) 



fe=i 



n3 



(1') ii') 

Taking into account Eqs. (176) and (177) (with L\-' = L^^' and of course also 



with a = and L^^^ = 0), Eq. (181) becomes a closed equation for L^^h 
'•^ fe=i fe=i 

(182) 

The physical root L^^) of Eq. (182) is the one vanishing in the HS limit Tjj — > 
00. Once known, Eq. (180) gives the contact values. 

This first-order approximation obtained from the RFA method turns out 
to coincide with the exact solution of the PY theory for SHS [95]. 



Second- Order Approximation 



As in the case of HS mixtures, a more flexible proposal is obtained by keeping 
a (and, conseq 
(178), one has 



a (and, consequently, L^ij^) different from zero. In that case, instead of Eq. 



27re'"«*Gy(s) = 



r(3) 

a 




+ 0(s-2). 



This implies 



(3) 



6tj 



7- (3) 



(183) 



(184) 



(185) 



If wc fix Vijiaij), Eqs. (176), (177), (184), and (185) allow one to express L(°), 
|_(i)^ L(2)^ a,nd L^^^ as linear functions of a. Thus, only the scalar parameter 
a remains to be fixed, analogously to what happens in the HS case. As done 
in the latter case, one possibility is to choose a in order to reproduce the 
isothermal compressibility x given by Eq. (148). To do so, one needs to find 
the coefficients H-j^ appearing in Eq. (152). The result is [104] 



H(0) = cW-(l-A(°))"\ 
H(i)=C«.(l-A(°))"\ 



(186) 
(187) 



where 
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N N N 

/s=l fc=l fe=l 



(188) 



N 

-E 

/c=l 



iV 



N 



N 



(3) 



.A 



fc=l 



i") (< 



kj 



kj 



E 

=1 

)■ 



fc=i 



l(0) 



o-(O) 



(189) 



(n + B)!^'^ 



(0) 



(1) 



n+1 



(nV2)!^'J 
(1) 



+ 



(2) 



(3) 



n! 



(190) 



Equation (187) gives H(i) in terms of a: hIj' = Py(a)/[Q(a)]2, where P^j(a) 
denotes a polynomial in a of degree 2A'^ and Q{a) denotes a polynomial of 
degree A^. It turns out then that, seen as a function of a, x is the ratio of 
two polynomials of degree 2N. Given a value of x, one may solve for a. The 
physical solution, which has to fulfill the requirement that Gij{s) is positive 
definite for positive real s, corresponds to the smallest positive real root. 

Once a is known, the scheme is complete: Eq. (184) gives L^^', then L^^^ is 
obtained from Eq. (185), and finally L^^^ and L^"^ arc given by Eqs. (176) and 
(177), respectively. Explicit knowledge of Gij{s) through Eqs. (154), (174), 
and (175) allows one to determine the Fourier transform hij{q) and the struc- 
ture factor Sij{q) through Eqs. (161) and (162), respectively. Finally, inverse 
Laplace transformation of Gij (s) yields Qij (r) [90] . 



Single Component SHS Fluids 

The special case of single component SHS fluids [105, 106] can be obtained 
from the multicomponent one by taking aij = a and Tij = r. Thus, the 
Laplace transform of rg{r) in the RFA is 

~ 27rs2 l + as-p [(^2(5)^^°) + Viis)L^^^ + Ms)L^^^ " e-«L(3)] ' 

(191) 

where we have taken a = 1. Equations (176) and (177) become 

The choice a = L^^^ = makes Eq. (191) coincide with the exact solution 
to the FY approximation for SHS [94], where L^^^ is the physical root (i.e., the 
one vanishing in the limit r ^ 00) of the quadratic equation [see Eq. (182)] 



48 M. Lopez de Haro, S. B. Yuste and A. Santos 



12tL(2) = 27r 




1^ 1 — n TT 



(194) 



We can go beyond the PY approximation by prescribing a contact value 
y(l), so that, according to Eqs. (184) and (185), 



By prescribing the isothermal compressibility %, the parameter a can be ob- 
tained as the physical solution (namely, the one remaining finite in the limit 
T oo) of a quadratic equation [106]. Thus, given an EOS for the SHS 
fluid, one can get the thermodynamically consistent values of y(l) and % and 
determine from them all the coefficients appearing in Eq. (191). 

Figure 12 shows the cavity function for -q = 0.164 and r = 0.13 as obtained 
from Monte Carlo simulations [101] and as predicted by the PY and RFA 
theories, the latter making use of the EOS recently proposed by Miller and 
Frenkel [102] . It can be observed that the RFA is not only more accurate than 
the PY approximation near r = 1 but also near r = 2. On the other hand, 
none of these two approximations account for the singularities (delta-peaks 
and/or discontinuities) of y{r) at r = ^/S/3, 5/3, Vs, 2, . . . [100, 101]. 



T 



(195) 




(196) 
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Fig. 12. Cavity function of a single component SHS fluid for 77 = 0.164 and r = 0.13. 

The solid line represents simulation data [101] . The dotted and dashed lines represent 
the PY and RFA approaches, respectively. 
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4.2 Single Component Square- Well Fluids 

Now we consider again the SW interaction potential (169) but for a single 
fluid, i.e., aij = a, Cij = e, Rij = R. Since no exact solution of the PY theory 
for the SW potential is known, the application of the RFA method is more 
challenging in this case than for HS and SHS fluids. 

As in the cases of HS and SHS, the key quantity is the Laplace transform 
of rg{r) defined by Eq. (98). It is again convenient to introduce the auxiliary 
function SI/{s) through Eq. (99). As before, the conditions g{r) = finite and 
X = finite imply Eqs. (102) and (104), respectively. However, the important 
difference between HS and SHS fluids is that in the latter case G{s) must 
reflect the fact that g{r) is discontinuous at r = i? as a consequence of the 
discontinuity of the potential ^(r) and the continuity of the cavity function 
y(r). This implies that G{s)^ and hence !f'(s), must contain the exponential 
term e~(^~'^)*. This manifests itself in the low-density limit, where the con- 
dition lim^^o y{i') = 1 yields 

1 

hm<lr{s) = ^gi/T.(i + _ g-(iJ-l)«(gl/T- _ l)(l + Rs) ' ^^^^^ 

where T* = ksT/e and we have taken a = 1. 

In the spirit of the RFA method, the simplest form that complies with Eq. 
(102) and is consistent with Eq. (197) is [107] 

2^1 + Qo + Oi5-c-(«-i)«(go + Q2s)' ^ ' 

where the coefficients Qo, Qi, Q2, Ei, E2, and E^ are functions of 77, T*, and 
R. The condition (104) allows one to express the parameters Qi, Ei, E2, and 
£^3 as linear functions of Qo and Q2 [107, 108] : 

Qi = fi + ? + MR' - i)Q2 - ^{R - if{R' + 2R + 3)Qo 

+Q2 - (i? - l)Qo, (199) 
Ex = [3 - 4(i?=^ - 1)Q2 + (i? - \f{R^ + 2ii + 3)Qo] , (200) 

^2 = T-Sr {1 - »? - 2(i? - 1) [1 - 2r]R{R + 1)] Q2 

+(i?-l)2[(l-r7(i?+l)2]go}, (201) 

= { (1 - + MR ~i){R+l- 2vR') Q2 

-T]{R - 1)2[4 + 2R- ri{3R^ + 2R+ l)]Qo} • (202) 
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Prom Eq. (102), we have 

5(1+) = §i. (203) 

The complete RDF is given by Eq. (122), where now Eq. (198) must be used 
in Eq. (123). In particular, tpi{r) and ilJ2{r) are 

Mr) = M{r)0{r) + ^Jnir + 1 - R)e{r + 1-R), (204) 

Mr) = M{r)O{r)+jP2i{r+l-R)O{r + l-R)+ij22{r+2-2R)0{r+2-2R), 

(205) 

where 

3 

V'u(r) = 2^E^|^«^^"^ (206) 

3 

M{r) = -477^X1 

Here, Sj are the three distinct roots of E{s) = —12r] + E^s + E^s^ + E^s^ and 

Ww{s) = l + Qo + Qxs, Wu{s) = -{Qo + Q2s). (208) 

W2o{s) = s[Ww{s)f, W2i{s) = 2sWw{s)Wii{s), W22{s) = s[Wii{s)]\ 

(209) 

To close the proposal, we need to determine the parameters Qq and Q2 by 
imposing two new conditions. An obvious condition is the continuity of the 
cavity function at r = i?, what implies 

g{R+) = e''^'g{R-). (210) 

This yields 

(1 _ e-i/r*^ ^^^(^ _ 1) = -^ii(o) = 27r^. (211) 

As an extra condition, we could enforce the continuity of the first derivative 
y' {r) iit r = R [109]. However, this complicates the problem too much without 
any relevant gain in accuracy. In principle, it might be possible to impose 
consistency with a given EOS, via either the virial route, the compressibility 
route, or the energy route. But this is not practical since no simple EOS for 
SW fluids is at our disposal for wide values of density, temperature, and range. 
As a compromise between simplicity and accuracy, we fix the parameter Qo 
at its exact zero-density limit value, namely Qo = e^/-^ — 1 [107]. Therefore, 
Eq. (211) becomes a transcendental equation for Q2 that needs to be solved 
numerically. For narrow SW potentials, however, it is possible to replace the 
exact condition (210) by a simpler one allowing Q2 to be obtained analytically 
[108] , which is especially useful for determining the thermodynamic properties 
[108,110]. 



rW2k{Si) + W^kis^)-W2k{s^) 



E'{s,) 



[E'{s,)Y 



(207) 
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Fig. 13. Radial distribution function of a single component SW fluid for R = 1.05, 
pa^ = 0.8, and T* = 0.5 (top panel), for R = 1.5, pa^ = 0.4, and T* = 1.5 (middle 
panel), and for R — 2.0, pa'* — 0.4, and T* = 3.0 (bottom panel). The circles 
represent simulation data [111] and the solid lines refer to the results obtained from 
the RFA method. 



It can be proven that the RFA proposal (198) reduces to the exact solutions 
of the PY equation [29,94] in the HS limit, i.e., e ^ or i? ^ 1, and in the 
SHS limit, i.e., e ^ oo and i? ^ 1 with {R - l)e^/'^' = finite [107, 108]. 

Comparison with computer simulations [107, 108, 110, 111] shows that the 
RFA for SW fluids is rather accurate at any fluid density if the potential well 
is sufficiently narrow (say R < 1.2), as well as for any width if the density 
is small enough (say p<T^ < 0.4). However, as the width and/or the density 
increase, the RFA predictions worsen, especially at low temperatures. As an 
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illustration, Fig. 13 compares the RDF provided by the RFA with Monte Carlo 
data [111] for three representative cases. 

4.3 Hard Disks 

As is well known, the PY theory is exactly solvable for HS fluids with an 
odd number of dimensions [112 114]. In particular, in the case of hard rods 
{d = 1), the PY theory provides the exact RDF g{r) or, equivalently, the exact 
cavity function y{r) outside the hard core (i.e., for r > a). However, it does 
not reproduce the exact y{r) in the overlapping region (i.e., for r < a) [85]. 
The full exact one- dimensional cavity function is [85] 

-(r-l)r,/(l-r)) °o ri''-^e-^r-n)ri/{l-v) 

,„,(.|„ = + g '(1 -,)-(„ -1). c - »)"*'(^ - »>. 

(212) 

where the subscript HR stands for hard rods and, as usual, a = 1 has been 
taken. Consequently, one has 

1 19 1 

5hr(1+I^) = I drr/iHR(r|r?) = 4%) = -- + -7?--y?^ (213) 

When d is even, the PY equation is not analytically solvable for the HS 
interaction. In particular, in the important case of hard disks (d = 2), one 
must resort to numerical solutions of the PY equation [1, 115]. Alternatively, 
a simple heuristic approach has proven to yield reasonably good results [116]. 
Such an approach is based on the naive assumption that the structure and 
spatial correlations of a hard-disk fluid share some features with those of a 
hard-rod and a hard-sphere fluid. This fuzzy idea becomes a more speciflc one 
by means of the following simple model [116]: 

5HD(r-|r?) = i/{r])guR{r\uJi{r])r]) + [1 - i^{v)]ms{r\i^3{v)v)- (214) 

Here, the subscript HD stands for hard disks {d = 2) and the subscript HS 
stands for hard spheres (rf = 3). The parameter 1^(77) is a density-dependent 
mixing parameter, while LOi{ri)ri and L03,{vi)r] are the packing fractions in one 
and three dimensions, respectively, which are "equivalent" to the packing 
fraction 77 in two dimensions. In Eq. (214), it is natural to take for gmiiA'n) 
the exact solution, Eq. (212). As for gn'R.{i'\'n)-: oiie might use the RFA recipe 
described in Section 3. However, in order to keep the model (214) as simple 
as possible, it is sufficient for practical purposes to take the PY solution, Eq. 
(108). In the latter approximation, 



, l + ??/2 p . , . „(0), ^ 10-277 + 77 
gHs(l+|?7) = n -77)2 ' drr/iHs(r-|7?) =ffiis'(??) = 



2 



{l-vr Jo 20(1 + 27?) 

(215) 
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In order to close the model (214), we still need to determine the parameters 
v{ri), LLii{rj), and ^3(77). To that end, wc first impose the condition that Eq. 
(214) must be consistent with a prescribed contact value gHD(l'''|»?) or, equiv- 
alently, with a prescribed compressibility factor Zhd(??) = 1 + 2r]gHr>{^~^\'il)^ 
with independence of the choice of the mixing parameter z^(r?). In other words, 

5HD(l+|r?) = 5hr(1+|wi(7?)??) = ffHs(l+|w3 (??)??)• (216) 
Making use of Eqs. (213) and (215), this yields 

gHD(l+|r/)-l 4ffHD(l+|r;) + 1 - V24gHD(l+|r;) + 1 

"^^^^= V9M1^\V) ' = 4,,HD(l+|r?) • 

(217) 

Once tOi{ri) and ^3(77) are known, we can determine 1/(7]) by imposing that 
the model (214) reproduces the isothermal compressibility Xhd('7) thermody- 
namically consistent with the prescribed Ziir,{r]) [cf. Eq. (97)]. From Eqs. (96) 
and (214) one has 

POO 

Xuuiv) = 1 + 877 / drr{v{r])hiin{r\uji{r])r]) + [1 - u{r])] husir\u^3iv)v)} , 
Jo 

(218) 

so that 

. . = {xMv)-l]/Sv-H^i{u;Mv) 

whore i7^]](r7) and i?Hs(??) are given by Eqs. (213) and (215), respectively. 

Once a sensible EOS for hard disks is chosen [see, for instance. Table 1], 
Eqs. (217) and (219) provide the parameters of the model (214). The results 
show that the scaling factor uJi{ri) is a decreasing function, while 1^3(7?) is 
an increasing function [116]. As for the mixing parameter !^(r?), it is hardly 
dependent of density and takes values around 1^(77) ~ 0.35-0.40. 

Comparison of the interpolation model (214) with computer simulation re- 
sults shows a surprisingly good agreement, despite the crudeness of the model 
and the absence of empirical fitting parameters, especially at low and mod- 
erate densities [116]. The discrepancies become important only for distances 
beyond the location of the second peak and for densities close to the stability 
threshold. 



5 Perturbation Theory 

When one wants to deal with realistic intcrmolccular interactions, the prob- 
lem of deriving the thermodynamic and structural properties of the system 
becomes rather formidable. Thus, perturbation theories of liquids have been 
devised since the mid twentieth century. In the case of single component flu- 
ids, the use of an accurate and well characterized RDF for the HS fluid in a 
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perturbation theory opens up the possibihty of deriving a closed theoretical 
scheme for the determination of the thermodynamic and structural proper- 
ties of more realistic models, such as the Lennard-Jones (LJ) fluid. In this 
section, we will consider this model system, which captures the basic physical 
properties of real non-polar fluids, to illustrate the procedure. 

In the application of the perturbation theory of liquids, the stepping stone 
has been the use of the HS RDF obtained from the solution to the PY equa- 
tion. Unfortunately, the absence of thermodynamic consistency present in the 
PY approximation (as well as in other integral equation theories) may clearly 
contaminate the results derived from its use within a perturbative treatment. 
In what follows we will reanalyze the different theoretical schemes for the ther- 
modynamics of L J fluids that have been constructed with perturbation theory, 
taking as the reference system the HS fluid. This includes the consideration 
of the RDF as obtained with the RFA method, which embodies thermody- 
namic consistency, as well as the proposal of a unifying framework in which 
all schemes fit in. With our development, we will be able to present a formula- 
tion which lends itself to relatively easy numerical calculations while retaining 
the merits that analytical results provide, namely a detailed knowledge and 
control of all the approximations involved. 

Let lis consider a three-dimensional fluid system defined by a pair inter- 
action potential (/'(r). The virial and energy EOS express the compressibility 
factor Z and the excess part of the Helmholtz free energy per unit volume 
f^^, respectively, in terms of the RDF of the system as 

Z=l-|7rp/3 rdr^giry, (220) 



r 



27rp/3 / dr<?^)(r)g(r)r^ (221) 
Jo 



pksT 

where /3 = l/fc^T. Let us now assume that 0(r) is split into a known (ref- 
erence) part (po{f) and a perturbation part 0i(r). The usual perturbative 
expansion for the Helmholtz free energy to first order in /? leads to [2] 



/ fo 



poo 

27rp/3 / drMr)9o{ry + O {(3^) , (222) 
Jo 



where /o and go{r) are the free energy and the RDF of the reference system, 
respectively. 

The LJ potential is 

</'L,i(r) =4e(r-i2_^-6)^ (223) 

where e is the depth of the well and, for simplicity, we have taken the distance 
at which the potential vanishes as the length unit, i.e., 0Lj(r = 1) = 0. For 
this potential the reference system may be forced to be a HS system, i.e., one 
can set 
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00 (r) = (/'Hs(r) 



00, r < (To, 
0, r > £70, 



(224) 



where ctq is a conveniently chosen effective HS diameter. In this case the 
Helmholtz free energy to this order is approximated by 



/] 



LJ 



/hs 



+ 2-Kpfi / dr(^Lj(r-)s'Hs(r'/o-o)r^ 



(225) 



Note that Eq. (225) may be rewritten in terms of the Laplace transform G{s) 
of (r/£7o)s'Hs(r-/o-o) as 



/] 



LJ 



/] 



HS 



where ^lj(s) satisfies 



pkeT 



2Trpl34 / ds^'Lj(s)G(s), 



{r)=ao dse-'-^/^°<l'Lj(s), 
Jo 



so that 



#Lj(s) = 4e(7o 



10! 



4! 



(226) 



(227) 



(228) 



Irrespective of the value of the diameter ao of the reference system, the 
right hand side of Eq. (226) represents always an upper bound for the value of 
the free energy of the real system. Therefore, it is natural to determine tro so 
as to provide the least upper bound. This is precisely the variational scheme of 
Mansoori and Canfield [117,118] and Rasaiah and Stell [119], usually referred 
to as MC/RS, and originally implemented with the PY theory for G{s), Eq. 
(108). In our case, however, we will consider G{s) as given by the RFA method, 
Eq. (113). Therefore, at fixed p and (3, the effective diameter ctq in the MC/RS 
scheme is obtained from the conditions 



d 

dao 



1 



+ 



10! 



10 



4! 



= 0, 



(229) 



V 



+ 



ds 0(51770) 



10! 



4! 



> 0. (230) 



In these equations, use has been made of the thermodynamic relationship 
between the free energy and the compressibility factor, Eq. (78). Moreover, we 
have called rjo = (7r/6)/3(jQ and have made explicit with the notation G(s|?7o) 
the fact that the HS RDF depends on the packing fraction rjQ. 
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Even if the reference system is not forced to be a HS fluid, one can still use 
Eq. (226) provided an adequate choice for (Tq is niadc such that the expansion 
involved in the right hand side of Eq. (222) yields the right hand side of Eq. 
(226) to order 0^. This is the idea of the Barker and Henderson [120] first 
order perturbation scheme (BHi), where the effective HS diameter is 



f 



ao= dr 



I _ g-/30Lj(r) 



(231) 



The same ideas may be carried out to higher order in the perturbation 
expansion. The inclusion of the second order term in the expansion yields the 
so-called macroscopic compressibility approximation [2] for the free energy, 
namely 



= ^^MrMry 



i 



7rp/3\o / dr <f>j {r)go {r^ + O {p') , (232) 



where xo is the (reduced) isothermal compressibility of the reference system 
[121]. 

To implement a particular perturbation scheme in this approximation un- 
der a unifying framework that eventually leads to easy numerical evaluation, 
two further assumptions may prove convenient. First, the perturbation poten- 
tial ^i(r) = (t)Lj{r) — 4'o{t) may be split into two parts using some "molecular 
size" parameter ^ > ctq such that 

Mr) = {tf}''z:-^' (233) 

Next, a choice for the RDF for the reference system is done in the form 

<?o(r)«0(r)yHs(r/ao), (234) 

where yus is the cavity (background) correlation function of the HS system 
and 6{r) is a step function defined by 

in which the functions 0a{r) and Ob{r) depend on the scheme. 

With these assumptions the integrals involved in Eq. (232) may be rewrit- 
ten as 

POO 

/„= / drcj>^{r)go{ry 
Jo 

r-cro 

= dr(t>Ur)Oa{r)yHsir/aoy+ drct>Ur)Oa{r)gHs{r/<Joy 

Jo Jan 



+ dr0?,(r)0fa(r)<?Hs(r/ao)r2, (236) 
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with n = 1,2 and where the fact that j/hsC^/c^o) = 5hs(''/o'o) when r > (Tq 
has been used. Decomposing the last integral as /^°° = — J^^ and applying 
the same step as in Eq. (226), Eq. (236) becomes 

roc (•(To 

Jo Jo 
+ dr [(f>Ur)Oa{r) - mOb{r)] 5Hs(r/ao)r^ (237) 

J era 

where the functions ^ife(s) and ^26 (s) are defined by the relation 

POO 

rcf>^f,{r)eb{r) = ao dse-^'^''°^„b{s). (238) 
Jo 

In the Barker-Henderson second order perturbation scheme (BH2), one 
takes 

Oa{r)=0, Ob{r) = l, ^ = ao, Mr)=0, ^ib{r) = {r~'^ - r-^) , 

(239) 

and (70 is computed according to Eq. (231). This choice ensures that 

= + 27: pP I dr (/)i(r)5Hs(r/o-o)r^ 

-TTpp\m / dr <l>l{r)gns{r/aoy + O {p^) . (240) 

J TO 

On the other hand, if one chooses 

0„(r)=exp[-/3(<^Lj(r)+e)], ^^(r) = 1, ^ = 2^/^, (241) 

<^i„(r) = -e, .^i6(r) = 4e (r-12 - r-*^) , (242) 

the scheme leads to the Weeks Chandler Andersen (WCA) theory [122] if one 
determines the HS diameter through the condition xo = Xhs [123], which in 
turn implies 



i 



I drr^e-^'^°('''yHs(r/ao) = / drrVs(r-/CTo) l-e'^'^^'M . (243) 

To close the scheme, the HS cavity function has to be provided in the range < 
r < ctq. Fortunately, relatively simple expressions for ynsif/^a) arc available 
in the literature [124-126], apart from our own proposal, Eq. (140). 

Note that 6'h(r) and (^i6(r), and thus also <l>nbis), are the same functions in 
the BH2 and WCA schemes. It is convenient, in order to have all the quantities 
needed to evaluate /lj in these schemes, to provide explicit expressions for 
<?ib(s) and (l>2b{s)- These are given by [cf. Eq. (228)] 



1/6 
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^16(s) =<5lj(.s), 



(244) 



101 



^26 (s) = 16eVo 



22! 16! 10! 



(245) 



Up to this point, we have embodied the most popular perturbation schemes 

within a unified framework that requires as input only the EOS of the HS fluid 
in order to compute the Helmholtz free energy of the LJ system and leads 
to relatively easy numerical computations. It should be clear that a variety 
of other possible schemes, requiring the same httle input, fit in our unified 
framework, which is based on the RFA method for (?hs(?'/<''o) and G(s). Once 
/lj has been determined, the compressibility factor of the LJ fluid at a given 
order of the perturbation expansion readily follows from Eqs. (222) or (232) 
through the thermodynamic relation 



Taking into account that the HS fluid presents a fluid-solid transition at 
a freezing packing fraction r/f ~ 0.494 [127] and a solid-fluid transition at 
a melting packing fraction ?7i„ ~ 0.54 [127], the fluid-solid and solid-fluid 
coexistence lines for the L J system may be computed from the values (p, T) 
determined from the conditions {-k /Q)pa-Q{p, T) = rjf and {-k /&)p(Jq{p, T) = rim, 
respectively, with the effective diameter <tq{p,T) obtained using any of the 
perturbative schemes. Similarly, admitting that there is a glass transition in 
the HS fluid at the packing fraction r/g ~ 0.56 [128], one can now determine 
the location of the liquid-glass transition line for the LJ fluid in the (p, T) 
plane from the simple relationship {-n /(S)puq{p,T) = rig. With a proper choice 
for .^HS) it has been shown [76,129,130] that the critical point, the structure, 
and the phase diagram (including a glass transition) of the LJ fluid may be 
adequately described with this approach. 



6 Perspectives 

In this chapter we have given a self-contained account of a simple (mostly 
analytical) framework for the study of the thermodynamic and structural 
properties of hard-core systems. Whenever possible, the developments have 
attempted to cater for mixtures with an arbitrary number of components (in- 
cluding polydisperse systems) and arbitrary dimensionality. We started con- 
sidering the contact values of the RDF because they enter directly into the 
EOS and arc required as input in the RFA method to compute the structural 
properties. With the aid of consistency conditions, we were able to devise var- 
ious approximate proposals which, when used in conjunction with a sensible 
choice for the contact value of the RDF of the single component fluid (re- 
quired in the formulation but otherwise chosen at will), have been shown to 




(246) 
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be in reasonably good agreement with simulation results and lead to accurate 
EOS both for additive and non-additive mixtures. Some aspects of the results 
that follow from the use of these EOS were illustrated by looking at demixing 
problems in these mixtures, including the far from intuitive case of a binary 
mixture of non-additive hard spheres in infinite dimensionality. 

After that, restricting ourselves to three-dimensional systems, we described 
the RFA method as applied to a single component hard-sphere fluid and to 
a multicomponcnt mixture of HS. Using this approach, we have been able to 
obtain explicit analytical results for the RDF, the direct correlation function, 
the static structure factor, and the bridge function, in the end requiring as 
input only the contact value of the RDF of the single component HS fluid 
(or equivalently its compressibility factor). One of the nice assets of the RFA 
approach is that it eliminates the thermodynamic consistency problem which 
is present in most of the integral equation formulations for the computation 
of structural quantities. Once again, when a sensible choice for the single 
component EOS is made, we have shown, through the comparison between 
the results of the RFA approach and simulation data for some illustrative 
cases, the very good performance of our development. Also, the use of the 
RFA approach in connection with some other related systems (sticky hard 
spheres, square-well fluids, and hard disks) has been addressed. 

The final part of the chapter concerns the use of HS results for more 
realistic intermolecular potentials in the perturbation theory of liquids. In 
this instance we have been able to provide a unifying scheme in which the 
most popular perturbation theory formulations may be expressed and which 
was devised to allow for easy computations. We illustrated this for a LJ fluid 
but it should be clear that a similar approach might be followed for other fluids 
and in fact it has recently been done in connection with the glass transition 
of hard-core Yukawa fluids [131]. 

Finally, it should be clear that there are many facets of the equilibrium and 
structural properties of hard-core systems that may be studied with a simi- 
lar approach but that up to now have not been considered. For instance, the 
generalizations of the RFA approach for systems such as hard hypersphcres, 
non-additive hard spheres, square- well mixtures, penetrable spheres [132], or 
the Jagla potential [133] appear as interesting challenges. Similarly, the ex- 
tension of the perturbation theory scheme to the case of LJ mixtures seems 
a worthwhile task. We hope to address some of these problems in the future 
and would be very much rewarded if some others were taken up by researchers 
who might find these developments also a valuable tool for their work. 

References 

1. J. A. Barker and D. Henderson, Rev. Mod. Phys. 48, 587 (1976). 

2. D. A. McQuarrie, Statistical Mechanics (Harper & Row, N. Y., 1976). 

3. H. L. Friedman, A Course in Statistical Mechanics (Prentice Hall, Englewood 
Cliffs, 1985). 



60 



M. Lopez de Haro, S. B. Yuste and A. Santos 



4. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, (Academic Press, 
London, 1986). 

5. J. L. Lebowitz and D. Zomick, J. Chcm. Phys. 54, 3335 (1971). 

6. J. T. Jenkins and F. Mancini, J. Appl. Mcch. 54, 27 (1987). 

7. C. Barrio and J. R. Solana, J. Chem. Phys. 115, 7123 (2001); 117, 2451(E) 
(2002). 

8. J. L. Lebowitz, Phys. Rev. A 133, 895 (1964). 

9. H. Reiss, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 31, 369 (1959); E. 
Helfand, H. L. Frisch, and J. L. Lebowitz, J. Chem. Phys. 34, 1037 (1961); J. 
L. Lebowitz, E. Helfand, and E. Praestgaard, J. Chem. Phys. 43, 774 (1965). 

10. M. J. Mandell and H. Reiss, J. Stat. Phys. 13, 113 (1975). 

11. Y. Rosenfeld, J. Chem. Phys. 89, 4272 (1988). 

12. M. Heying and D. S. Corti, J. Phys. Chem. B 108, 19756 (2004). 

13. T. Boublik, J. Chem. Phys. 53, 471 (1970). 

14. E. W. Grundke and D. Henderson, Mol. Phys. 24, 269 (1972). 

15. L. L. Lee and D. Levesque, Mol. Phys. 26, 1351 (1973). 

16. G. A. Mansoori, N. F. Carnahan, K. E. Starling, and J. T. W. Leland, J. Chem. 
Phys. 54, 1523 (1971). 

17. D. Henderson, A. Malijevsky, S. Labik, and K. Y. Chan, Mol. Phys. 87, 273 

(1996). 

18. D. H. L. Yau, K.-Y. Chan, and D. Henderson, Mol. Phys. 88, 1237 (1996); 91, 
1137 (1997). 

19. D. Henderson and K. Y. Chan, J. Chem. Phys. 108, 9946 (1998); Mol. Phys. 
94, 253 (1998); 98, 1005 (2000). 

20. D. Henderson, D. Boda, K. Y. Chan, and D. T. Wasan, Mol. Phys. 95, 131 

(1998). 

21. D. Matyushov, D. Henderson, and K.-Y. Chan, Mol. Phys. 96, 1813 (1999). 

22. D. Cao, K.-Y. Chan, D. Henderson, and W. Wang, Mol. Phys. 98, 619 (2000). 

23. D. V. Matyushov and B. M. Ladanyi, J. Chem. Phys. 107, 5815 (1997). 

24. C. Barrio and J. R. Solana, J. Chem. Phys. 113, 10180 (2000). 

25. D. Viduna and W. R. Smith, Mol. Phys. 100, 2903 (2002); J. Chem. Phys. 
117, 1214 (2002). 

26. D. Henderson, Mol. Phys. 30, 971 (1975). 

27. A. Santos, M. Lopez de Haro, and S. B. Yuste, J. Chem. Phys. 103, 4622 
(1995); M. Lopez de Haro, A. Santos, and S. B. Yuste, Eur. J. Phys. 19, 281 

(1998) . 

28. S. Luding, Phys. Rev. E 63, 042201 (2001); S. Luding, Adv. Compl. Syst. 

4, 379 (2002); S. Luding and O. StrauB, in Granular Gases, T. Poschel and 

5. Luding, eds. (LNP 564, Springer- Verlag, Berlin, 2001), pp. 389-409. 

29. M. S. Wertheim, Phys. Rev. Lett. 10, 321 (1963); E. Thiele, J. Chem. Phys. 
39, 474 (1963). 

30. N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969). 

31. M. Luban and J. P. J. Michels, Phys. Rev. A 41, 6796 (1990). 

32. E. Hamad, J. Chem. Phys. 101, 10195 (1994). 

33. C. Vega, ,J. Chem. Phys. 108, 3074 (1998). 

34. N. M. Tukur, E. Z. Hamad, and G. A. Mansoori, J. Chem. Phys. 110, 3463 

(1999) . 

35. A. Santos, S. B. Yuste, and M. Lopez de Haro, Mol. Phys. 96, 1 (1999). 

36. A. Malijevsky and J. Veverka, Phys. Chem. Chem. Phys. 1, 4267 (1999). 



Alternative Approaches to Hard-Sphere Liquids 



61 



37. A. Santos, S. B. Yuste, and M. Lopez de Haro, Mol. Phys. 99, 1959 (2001). 

38. M. Gonzalez-Mclchor, J. Alejandre, and M. Lopez de Haro, J. Chem. Phys. 
114, 4905 (2001). 

39. M. Lopez de Haro, S. B. Yuste, and A. Santos, Phys. Rev. E 66, 031202 
(2002). 

40. A. Santos, Mol. Phys. 96, 1185 (1999); 99, 617(E) (2001). 

41. C. Regnaut, A. Dyan, and S. Amokrane, Mol. Phys. 99, 2055 (2001); 100, 
2907(E) (2002). 

42. A. Santos, S. B. Yuste, and M. Lopez de Haro, J. Chem. Phys. 117, 5785 
(2002). 

43. A. Santos, S. B. Yuste and M. Lopez de Haro, J. Chem. Phys. 123, 234512 

(2005) ; M. Lopez de Haro, S. B. Yuste, and A. Santos, Mol. Phys. 104, 3461 

(2006) . 

44. S. Luding and A. Santos, J. Chem. Phys. 121, 8458 (2004). 

45. M. Barosova, A. Malijevsky, S. Labik, and W. R. Smith, Mol. Phys. 87, 423 

(1996). 

46. H. Hansen-Goos and R. Roth, J. Chem. Phys. 124, 154506 (2006). 

47. R. Evans, in Liquids and Interfaces, edited by J. Charvolin, J. F. Joanny, and 
J. Zinn- Justin (North-Holland, Amsterdam, 1990). 

48. Y. Roscnfcld, Phys. Rev. Lett. 63, 980 (1989). 

49. A. Malijevsky, M. Barosova, and W. R. Smith, Mol. Phys. 91, 65 (1997). 

50. Al. Malijevsky, A. Malijevsky, S. B. Yuste, A. Santos, and M. Lopez de Haro, 
Phys. Rev. E 66, 061203 (2002). 

51. M. Buzzacchi, I. Pagonabarraga, and N. B. Wilding, J. Chem. Phys. 121, 11362 

(2004) . 

52. Al. Malijevsky, S. B. Yuste, A. Santos, and M. Lopez de Haro, preprint arXiv: 
cond-mat/0702284. 

53. F. Lado, Phys. Rev. E 54, 4411 (1996). 

54. I. Prigogine and S. Lafleur, Bull. Classe Sci. Acad. Roy Belg. 40, 484, 497 
(1954). 

55. S. Asakura and F. Oosawa, J. Chem. Phys. 22, 1255 (1954); J. Polym. Sci. 33, 
183 (1958). 

56. R. Kikuchi, J. Chem. Phys. 23, 2327 (1955). 

57. P. Ballone, G. Pastore, G. GaUi, and D. Gazzillo, Mol. Phys. 59, 275 (1986). 

58. D. Gazzillo, G. Pastore, and S. Enzo, J. Phys.: Condens. Matter 1, 3469 (1989); 
D. Gazzillo, G. Pgistore, and R. Frattini, J. Phys.: Condens. Matter 2,8465 
(1990). 

59. A. Santos, M. Lopez de Haro, and S. B. Yuste, J. Chem. Phys. 122, 024514 

(2005) . 

60. E. Z. Hamad, J. Chem. Phys. 105, 3229 (1996). 

61. E. Z. Hamad, J. Chem. Phys. 105, 3222 (1996). 

62. H. Hammawa and E. Z. Hamad, J. Chem. Soc. Faraday Trans. 92, 4943 (1996). 

63. M. Al-Naafa, J. B. El-Yakubu, and E. Z. Hamad, Fluid Phase Equilibria 154, 

33 (1999). 

64. J. Jung, M. S. Jhon, and F. H. Ree, J. Chem. Phys. 100, 528 (1994). 

65. J. Jung, M. S. Jhon, and F. H. Ree, J. Chem. Phys. 100, 9064 (1994). 

66. T. Coussaert and M. Bans, J. Chem. Phys. 109, 6012 (1998). 

67. A. Yu. Vlasov and A. J. Masters, Fluid Phase Equilibria 212, 183 (2003). 

68. M. Lopez de Haro and C. F. Tejero, J. Chem. Phys. 121, 6918 (2004). 



62 M. Lopez de Haro, S. B. Yuste and A. Santos 



69. S. B. Yuste, A. Santos, and M. Lopez de Haro, Europhys. Lett. 52, 158 (2000). 

70. H.-O. Carmosin, H. L. Frisch, and J. K. Percus, J. Stat. Phys. 63, 791 (1991). 

71. A. Santos and M. Lopez de Haro, Phys. Rev. E 72, 010501(R) (2005). 

72. R. Roth, R. Evans, and A. A. Louis, Phys. Rev. E 64, 051202 (2001). 

73. S. B. Yuste and A. Santos, Phys. Rev. A 43, 5418 (1991). 

74. S. B. Yuste, M. Lopez de Haro, and A. Santos, Phys. Rev. E 53, 4820 (1996). 

75. M. Roblcs, M. Lopez de Haro, A. Santos, and S. B. Yuste, J. Chem. Phys. 108, 
1290 (1998). 

76. M. Robles and M. Lopez de Haro, Europhys. Lett. 62, 56 (2003). 

77. M. Robles and M. Lopez de Haro, J. Chem. Phys. 107, 4648 (1997). 

78. E. Waisman, Mol. Phys. 25, 45 (1973); D. Henderson and L. Blum, Mol. Phys. 
32, 1627 (1976); J. S. H0yc and L. Blum, J. Stat. Phys. 16, 399 (1977). 

79. A. Di'ez, J. Largo, and J. R. Solana, J. Chem. Phys. 125, 074509 (2006). 

80. J. Kolafa, S. Labik, and A. Malijevsky, Phys. Chem. Chem. Phys. 6, 2335 
(2004). See also http://www.vscht.cz/fch/software/hsmd/ for molecular dy- 
namics results of g{r). 

81. A. Trokhymchuk, I. Nezbeda, J. Jirsak, and D. Henderson, J. Chem. Phys. 
123, 024501 (2005). 

82. M. Lopez de Haro, A. Santos, and S. B. Yuste, J. Chem. Phys. 124, 236102 

(2006). 

83. L. L. Lee, J. Chem. Phys. 103, 9388 (1995); L. L. Lee, D. Ghonasgi, and E. 
Lomba, J. Chem. Phys. 104, 8058 (1996); L. L. Lee and A. Malijevsky, J. 
Chem. Phys. 114, 7109 (2001). 

84. S. Labfk and A. Malijevsky, Mol. Phys. 53, 381 (1984). 

85. Al. Malijevsky and A. Santos, J. Chem. Phys. 124, 074508 (2006). 

86. A. Santos and Al. Malijevsky, Phys. Rev. E 75, 021201 (2007). 

87. S. B. Yuste, A. Santos, and M. Lopez de Haro, J. Chem. Phys. 108, 3683 
(1998). 

88. L. Blum and J. S. H0ye, J. Phys. Chem. 81, 1311 (1977). 

89. J. Abate and W. Whitt, Queuing Systems 10, 5 (1992). 

90. A code using the Mathematica computer algebra system to obtain Gij{s) 
and gij{r) with the present method is available from the web page 
http://www.unex.es/eweb/fisteor/santos/filesRFA.html. 

91. N. W. Ashcroft and D. C. Langreth, Phys. Rev. 156, 685 (1967). 

92. A. B. Bathia and D. E. Thornton, Phys. Rev. B 8, 3004 (1970). 

93. S. B. Yuste, A. Santos, and M. Lopez de Haro, Mol. Phys. 98, 439 (2000). 

94. R. J. Baxter, J. Chem. Phys. 49, 2270 (1968). 

95. ,J. W. Pcrram and E. R. Smith, Chem. Phys. Lett. 35, 138 (1975). 

96. B. Barboy Chem. Phys. 11, 357 (1975); B. Barboy and R. Tcnnc, Chem. Phys. 
38, 369 (1979). 

97. G. Stell, J. Stat. Phys. 63, 1203 (1991); B. Borstnik, C. G. Jesudason, and G. 
Stell, J. Chem. Phys. 106, 9762 (1997). 

98. B. Barboy J. Chem. Phys. 61, 3194 (1974). 

99. J. W. Perram and E. R. Smith, Chem. Phys. Lett. 39, 328 (1975); P. T. 
Cummings, J. W. Pcrram, and E. R. Smith, Mol. Phys. 31, 535 (1976); E. R. 
Smith and J. W. Perram, J. Stat. Phys. 17, 47 (1977); J. W. Perram and E. 
R. Smith, Proc. R. Soc. London A353, 193 (1977); W. G. T. Kranendonk and 
D. Frenkel, Mol. Phys. 64, 403 (1988); C. Regnaut and J. C. Ravey J. Chem. 
Phys. 91, 1211 (1989); G. Stell and Y. Zhou, J. Chem. Phys. 91, 3618 (1989); 



Alternative Approaches to Hard-Sphere Liquids 



63 



J. N. Herrera and L. Blum, J. Chem. Phys. 94, 6190 (1991); A. Jamnik, D. 
Bratko, and D. J. Henderson, J. Chem. Phys. 94, 8210 (1991); S. V. G. Mcnon, 
C. Manohar, and K. S. Rao, J. Chem. Phys.; 95, 9186 (1991); Y. Zhou and G. 
Stall, J. Chem. Phys. 96, 1504 (1992); E. Dickinson, J. Chom. Soc. Faraday 
Trans. 88, 3561 (1992); C. F. Tcjero and M. Bans, Phys. Rev. E, 48, 3793 
(1993); K. Shuklaand R. Rajagopalan, Mol. Phys. 81, 1093 (1994); C. Regnaut, 
S. Amokrane, and Y. Hcno, J. Chem. Phys. 102, 6230 (1995); C. Regnaut, S. 
Amokrane, and P. Bobola, Prog. Colloid Polym. Sci. 98, 151 (1995); Y. Zhou, 
C. K. Hall, and G. Stell, Mol. Phys. 86, 1485 (1995); J. N. Herrera^Pacheco 
and J. F. Rojas-Rodn'guez, Mol. Phys. 86, 837 (1995); Y. Hu, H. Liu, and J. M. 
Prausnitz, J. Chem. Phys. 104, 396 (1996); O. Bernard and L. Blum, J. Chem. 
Phys. 104, 4746 (1996); L. Blum, M. F. Holovko, and L A. Protsykevych, J. 
Stat. Phys. 84, 191 (1996); S. Amokrane, P. Bobola and C. Regnaut, Prog. 
Colloid Polym. Sci. 100, 186 (1996); S. Amokrane and C. Regnaut, J. Chem. 
Phys. 106, 376 (1997); C. Tutschka, G. Kahl, and E. Ricgler, Mol. Phys. 100, 
1025 (2002); D. Gazzillo and A. Giacometti, Mol. Phys. 100, 3307 (2002); M. 
A. Miller and D. Frenkel, Phys. Rev. Lett. 90, 135702 (2003); D. Gazzillo and 
A. Giacometti, J. Chem. Phys. 120, 4742 (2004); R. Fantoni, D. Gazzillo, and 
A. Giacometti, Phys. Rev. E 72, 011503 (2005); A. Jamnik, Chem. Phys. Lett. 
423, 23 (2006). 

100. A. J. Post and E. D. Glandt, J. Chem. Phys. 84, 4585 (1986); N. A. Seaton 
and E. D. Glandt, J. Chem. Phys. 84, 4595 (1986); 86, 4668 (1986); 87, 1785 
(1987). 

101. M. A. Miller and D. Frenkel, J. Phys.: Condens. Matter 16, S4901 (2004). 

102. M. A. Miller and D. Frenkel, J. Chem. Phys. 121, 535 (2004). 

103. Al. Malijevsky, S. B. Yuste, and A. Santos, J. Chem. Phys. 125, 074507 (2006). 

104. A. Santos, S. B. Yuste, and M. Lopez de Haro, J. Chem. Phys. 109, 6814 
(1998). 

105. S. B. Yuste and A. Santos, J. Stat. Phys. 72, 703 (1993). 

106. S. B. Yuste and A. Santos, Phys. Rev. E 48, 4599 (1993). 

107. S. B. Yuste and A. Santos, J. Chem. Phys. 101, 2355 (1994). 

108. L. Acedo and A. Santos, J. Chem. Phys. 115, 2805 (2001). 

109. L. Acedo, J. Stat. Phys. 99, 707 (2000). 

110. J. Largo, J. R. Solana, L. Acedo, and A. Santos, Mol. Phys. 101, 2981 (2003). 

111. J. Largo, J. R. Solana, S. B. Yuste, and A. Santos, J. Chem. Phys. 122, 084510 
(2005). 

112. C. Freasier and D. J. Isbister, Mol. Phys. 42, 927 (1981). 

113. E. Leuthcusscr, Physica A 127. 667 (1984). 

114. M. Roblos, M. Lopez dc Haro, and A. Santos, J. Chem. Phys. 120, 9113 (2004). 

115. D. G. Chae, F. H. Ree, and T. Ree, J. Chem. Phys. 50, 1581 (1976). 

116. S. B. Yuste and A. Santos, J. Chem. Phys. 99, 2020 (1993). 

117. G. A. Mansoori and F. B. Canfield, J. Chem. Phys. 51, 4958 (1969). 

118. G. A. Mansoori, J. A. Provine, and F. B. Canfield, J. Chem. Phys. 51, 5295 
(1969). 

119. J. Rasaiah and G. StcU, Mol. Phys. 18, 249 (1970). 

120. J. A. Barker and D. Henderson, J. Chem. Phys. 47, 2856 (1967). 

121. The macroscopic compressibility approach is only one of the possibilities of ap- 
proximation to the second order Barker-Henderson perturbation theory term. 
Another successful approach is the local-compressibility approximation (see 



64 M. Lopez de Haro, S. B. Yuste and A. Santos 



Ref. [2], p. 308). This expresses the free energy in terms of 0i(r') and HS 
quantities. 

122. J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chem. Phys. 53, 149 (1971). 

123. A simple algorithm to compute a rather accurate approximation for the HS 
diameter ao in the WCA theory has been given in L. Verlet and J. J. Weis, 
Phys. Rev. A 5, 939 (1972). 

124. D. Henderson and E. W. Grundkc, J. Chcm. Phys. 63, 601 (1975). 

125. J. A. Ballance and R. J. Speedy, Mol. Phys. 54, 1035 (1985). 

126. Y. Zhou and G. Stell, J. Stat. Phys. 52, 1389 (1988). 

127. J.-P. Hansen and L. Verlet, Phys. Rev. 184, 151 (1969). 

128. R. J. Speedy, J. Chem. Phys. 100, 6684 (1994). 

129. M. Robles and M. Lopez dc Haro, Phys. Chcm. Chcm. Phys. 3. 5528 (2001). 

130. M. Lopez de Haro and M. Robles, J. Phys.: Condens. Matt. 16, S2089 (2004). 

131. M. Lopez de Haro and M. Robles, Physica A 372, 307 (2006). 

132. C. N. Likos, Phys. Rep. 348, 267 (2001). 

133. E. A. Jagla, J. Chem. Phys. Ill, 8980 (1999). 



